This article is mentioned in:

Introduction

Skin cancer is a malignant tumor that afflicts individuals all over the world. Skin cancer is divided into malignant melanoma and non-melanoma skin cancer (NMSC) (1,2). According to global cancer statistics in 2018, there were 142,056 new cases of NMSC, accounting for 5.8% of global cancer cases, and 65,155 NMSC-related deaths, accounting for 0.7% of global cancer mortality. There were 779,723 new cases of skin malignant melanoma, accounting for 1.6% of global cancer cases, and 60,712 skin melanoma-related deaths, accounting for 0.6% of global cancer mortality (3).

According to the 2015 Chinese Cancer Statistics Report, the incidence of skin cancer in China was 8/1,000 and the mortality rate was 3.22/1,000 (4,5). Compared with other types of cancer, skin cancer has a higher recurrence rate, with a 35% probability of recurrence in the first 3 years and a 50% probability of recurrence in the first 5 years. Recurrent skin cancer is usually the same sub-type as the original cancer (6). Although the mortality rate of skin cancer is extremely low and the cure rate is high, its high incidence and recurrence rate constitute a significant economic burden to health services. Furthermore, skin cancer located on the head and other highly visible areas may affect a patients mental wellbeing and quality of life (7,8).

Cutaneous metastases of malignant tumors can be caused by malignant tumor cells traveling through the blood or lymphatic system, tissue interstitial diffusion (tissue gap diffusion), or surgical implantation (9). The higher the malignancy of the tumor, the more likely it is to metastasize, and skin metastasis is often the end-stage manifestation of a malignant tumor (10). Once metastasis occurs, the prognosis of cancer is poor (11). Unfortunately, skin metastasis may be the first clinical manifestation of a malignant tumor (12). Therefore, a suitable prognosis serves an important role in the recurrence and treatment of metastatic skin cancer. The aim of the present study was to identify potential prognostic biomarkers of metastatic skin cancer, which may be used in a clinical setting, through data mining and analysis of skin cancer prognosis genes.

With the development of gene chip technology and next-generation sequencing technology (13), and the constant revision of the viewpoints of individualized medical treatment and precision medical care, understanding skin cancer from the perspective of the genome and proposing more effective genetic biomarkers provides more relevant information for drug development and clinical decision-making (14-18). Traditional statistical methods often yield unstable results and excessive errors when applied to gene expression analysis (19). Furthermore, in high-throughput gene expression experiments, the number of variables is generally much higher than the sample size, which is called the Curse of Dimensionality (20). Therefore, in the present study, least absolute shrinkage and selection operator (LASSO) was used to mine genomic data to minimize the instability caused by high-dimensional data (21-27).

Materials and methods

Data collection

Using The Cancer Genome Atlas, data on skin cancer from the Xena Functional Genomics Explorer (xenabrowser.net/datapages/) was obtained, including the number of patients (n=481) and the number of genes assessed (n=20,530). Among these, the number of primary samples (Primary Tumor) was 105, the number of metastatic (Metastatic) samples was 368, and the number of other types of samples was 8 (28-32).

Data preprocessing

First, 219 samples with missing survival attribute values were removed from 473 primary and metastatic samples and 255 patient samples were retained. These data were stratified into primary tumor samples (n=72) and metastatic samples (n=183) by matching the patients' samples in the gene expression spectrum. Finally, using the ‘sample’ function in R version 3.5.2(33), the metastatic samples were randomly divided into two groups; training samples (n=91) and test samples (n=92).

Clinical data analysis

Several clinical factors affect the prognosis of patients with skin cancer, so it is necessary for data analysis to consider multiple attribute values in the sample as influencing factors. Based on previous studies, clinical variables including age, sex, history of radiotherapy, Tumor-Node-Metastasis pathological stage (34) and cancer status were used to analyze the related clinical factors in the subsequent data mining analysis (35-38).

Analysis of gene data

First, the ‘Limma’ package version 3.42.2(39) in R was used to analyze the differentially expressed genes of the 91 training set samples and 72 primary tumor samples. According to the adjusted P-value (adj.P.val<0.001), 783 genes were considered to be differentially expressed genes.

The ‘survival’ package (rdocumentation.org/packages/survival; version 3.1-8) in R was used to perform Cox regression analysis of the differentially expressed genes and the Cox coefficient, the hazard ratio (HR) and the P-values of the Wald test of each gene were calculated using the Kaplan-Meier method (40), and the genes that were significantly associated with the survival of the patients (P<0.05) were screened using the ‘survivalROC’ package (version number 1.0.3; rdocumentation.org/packages/survivalROC).

Finally, the screened genes were analyzed again using LASSO in the R package ‘glmnet’ (rdocumentation.org/packages/glmnet; version 3.0-2), to obtain more critical genes. Through 10-fold cross validation, 20 risk genes, which were closely related to survival, were identified (41-49). Additionally, the Gene Ontology (GO) Cell Component Ontology Method (50,51) was used to analyze pathway involved in protein activity.

Prognostic index (PI) calculation

As an important indicator of the integration of risk genes, a PI value can be determined for each patient with skin cancer. The PI was obtained by linearly fitting the product of the expression and the coefficient corrected by LASSO of each gene. The formula for calculating the PI was: PI=β1 X1+β2X2+…+βiXi…+βnXn; where Xi is the expression of the ith gene and βi is the coefficient corrected by LASSO of the ith gene.

Data validation

By extrapolating genes with a P<0.05, the 20 risk genes obtained were verified using the test sample using the same methods described above.

Statistical analysis of clinical variables

From the clinical information of 255 patient samples, the patients' age, sex, radiation therapy, pathological T-stage and cancer status were taken as single variables. Cox regression analysis was performed, and the P-value of the Log-rank test and the HR value of each clinical factor were calculated. P<0.05 was considered to indicate a statistically significant difference.

Results

Clinicopathological variables

Age, pathologic T-stage and cancer status were all significantly associated with the recurrence of skin cancer, suggesting that these three clinical factors may be used as independent prognostic indicators (Table I).

Table I Clinicopathological characteristics of the patients.

Table I

Clinicopathological characteristics of the patients.

Characteristics	Number of patients (relapse)	P-value	HR (95% CI)
Age		0.030	1.669 (1.037-2.688)
>57	123(35)
≤57	132(36)
Sex		0.200	1.359 (0.8058-2.291)
Male	160(49)
Female	95(22)
Cancer status
Unknown	3(1)	7x10-10
With tumor	160(18)		0.802 (0.1089-5.908)
Tumor free	92(52)		0.1572 (0.0206-1.201)
Radiation therapy		0.200	1.7204 (0.851-3.478)
Unknown	1 (0)
Yes	38(9)
No	216(62)
Pathological T-stage		0.009	2.003 (1.175-3.412)
T1-T3	165(49)
T4	90(22)

[i] HR, hazard ratio; CI, confidence interval.

Genetic data analysis

After determining differential gene expression, 44 genes were considered significantly associated with the survival of patients using survival analysis. Cox regression analysis and LASSO analysis were performed, and 20 risk genes associated with skin cancer were identified (Table II).

Table II Prognostic genes.

Table II

Prognostic genes.

Gene symbol	Name	Univariate HR (95% CI)	Coefficient	P-value	LASSO HR (95% CI)
TMEM45B	Transmembrane protein 45B	0.899 (0.822-0.982)	-0.086	1.88x10-2	0.918 (0.864-0.975)
CDKN1B	Cyclin-dependent kinase inhibitor 1B (p27, Kip1)	0.928 (0.886-0.972)	-0.029	1.55x10-3	0.971 (0.914-1.032)
PKHD1L1	Polycystic kidney and hepatic disease 1 (autosomal recessive)-like 1	0.938 (0.886-0.993)	-0.010	2.73x10-2	0.990 (0.932-1.052)
PLCL2	Phospholipase C-like 2	0.961 (0.929-0.994)	-0.006	2.14x10-2	0.994 (0.936-1.056)
CKMT1A	Creatine kinase, mitochondrial 1A	1.024 (1.000-1.049)	0.001	4.90x10-2	1.001 (0.942-1.063)
SPINK5	Serine peptidase inhibitor, Kazal type 5	1.027 (1.002-1.052)	0.007	3.31x10-2	1.008 (0.948-1.070)
CST6	Cystatin E/M	1.027 (1.000-1.054)	0.011	4.67x10-2	1.011 (0.952-1.074)
FABP5	Fatty acid binding protein 5 (psoriasis-associated)	1.037 (1.007-1.067)	0.007	1.53x10-2	1.007 (0.948-1.070)
PTK6	Protein tyrosine kinase 6	1.047 (1.010-1.084)	0.029	1.11x10-2	1.030 (0.969-1.094)
TXNDC17	Thioredoxin domain containing 17	1.047 (1.002-1.095)	0.004	4.02x10-2	1.004 (0.945-1.066)
LTB4R	Leukotriene B4 receptor	1.050 (1.007-1.096)	0.018	2.39x10-2	1.018 (0.958-1.082)
FAM100A	Family with sequence similarity 100, member A	1.050 (1.008-1.094)	0.015	1.90x10-2	1.016 (0.956-1.079)
KRT10	keratin 10	1.057 (1.025-1.091)	0.016	4.33x10-4	1.016 (0.956-1.080)
PPP1R13L	Protein phosphatase 1, regulatory subunit 13 like	1.059 (1.010-1.109)	0.000	1.65x10-2	1.000 (0.941-1.063)
DLX3	Distal-less homeobox 3	1.063 (1.025-1.102)	0.008	9.40x10-4	1.008 (0.949-1.071)
KPRP	Keratinocyte proline rich protein	1.064 (1.015-1.115)	0.044	9.68x10-3	1.045 (0.984-1.111)
VSIG10L	V-set and immunoglobulin domain containing 10 like	1.066 (1.011-1.123)	0.021	1.87x10-2	1.021 (0.961-1.085)
MYOM3	Myomesin 3	1.073 (1.015-1.134)	0.040	1.25x10-2	1.041 (0.979-1.106)
NMU	Neuromedin U receptor 1	1.092 (1.042-1.144)	0.072	2.22x10-4	1.075 (1.012-1.142)
PIGW	Phosphatidylinositol glycan anchor biosynthesis, class W	1.113 (1.044-1.188)	0.002	1.08x10-3	1.002 (0.943-1.064)

[i] HR, hazard ratio; CI, confidence interval.

Gene PI analysis

Through linear fitting of the product of expression and regression coefficient of the 20 genes in each sample, the PI of each patient was calculated, and the patients were sorted from lowest to highest according to their PI value. Based on the median PI value, the patients were divided into high-risk and low-risk groups (Fig. 1). The lower the PI value, the lower the risk of recurrent survival, and the higher the PI value, the higher the risk of recurrent survival.

Figure 1 Stratification of patients with skin cancer into high-risk and low-risk groups using the prognostic gene biomarkers. PI, prognostic index.

Using the Kaplan Meier method, the survival curves for the two groups of patients, are presented in Fig. 2A; patients considered low risk exhibited significantly longer overall survival times (P<0.001; HR=26.321). Using the 5-year survival rates, a Receiver Operating Characteristic (ROC) curve was drawn (Fig. 2B) and analysis was performed using the ‘survivalROC’ package. The advantages and disadvantages of the model constructed using the 20 gene biomarkers were determined based on the Area Under the Curve (AUC) value. The results showed that AUC was equal to 0.908 (AUC >0.5 indicates a suitable model). The heat map of risk gene expression profiles in high-risk and low-risk patients were plotted (Fig. 3). The high-risk group was clearly distinguished from the low-risk group. This indicated that the models constructed by these 20 gene biomarkers performed well.

Figure 2 Integrative model for predicting outcome. (A) Survival curves of the high-risk and low-risk patients distinguished according to the prognostic index value. There was a significant difference in survival between the two groups. (B) ROC curve analysis of gene marker models. ROC, Receiver Operating Characteristic; AUC, area under the curve.

Figure 3 Heat map of the expression profile of the risk genes in the high-risk and low-risk patients.

Based on the experimental results, the 20 gene prognostic biomarkers could be used to significantly classify the patients with skin cancer in the training samples into two groups: High risk and low risk. In order to further verify the accuracy of the test results, the 20 genetic biomarkers were used to validate the test samples. As shown in Fig. 4, these genetic biomarkers could still classify patients with skin cancer in the test sample into high-risk and low-risk categories (P=0.004, HR=1.194). The AUC of the ROC curve was 0.855.

Figure 4 Integrative model for validating the results. (A) Survival curves of patients in the validation set stratified by the PI value. (B) ROC curve analysis of gene marker models. ROC, Receiver Operating Characteristic; AUC, area under the curve.

Prognostic gene function analysis

The online genetic analysis tool STRING (52) was used to analyze and study the association between the identified gene biomarkers and their associated protein synthesis pathways (Fig. 5; Table III) Gene Ontology (GO) cellular component ontology analysis identified a pathway of which involved a cyclin-dependent protein kinase holoenzyme complex, which included three proteins (CCND1, CDK2 and CDK4).

Figure 5 Protein network of CDK4, CDKN18, CCNA2, CCND1, PTK6 and CDK2.

Table III GO biological processes associated with the prognostic genes.

Table III

GO biological processes associated with the prognostic genes.

Pathway ID	Pathway description	Gene count	False discovery rate	Matching proteins
GO.0010948	Negative regulation of cell cycle process	5	0.0219	CCNA2, CCND1, CDK2, CDK4, CDKN1B
GO.0031100	Organ regeneration	4	0.0219	CCNA2, CCND1, CDK2, CDK4
GO.0044773	mitotic DNA Damage checkpoint	4	0.0219	CCNA2, CCND1, CDK2, CDKN1B
GO.0071156	Regulation of cell cycle arrest	4	0.0219	CCND1, CDK2, CDK4, CDKN1B
GO.1901990	Regulation of mitotic cell cycle phase transition	5	0.0219	CCNA2, CCND1, CDK2, CDK4, CDKN1B
GO.0032355	Response to estradiol	4	0.0238	CCNA2, CCND1, CDK2, CDKN1B
GO.0010389	Regulation of G2/M transition of mitotic cell cycle	3	0.0277	CCNA2, CCND1, CDK4
GO.0044772	Mitotic cell cycle phase transition	5	0.0277	CCNA2, CCND1, CDK2, CDK4, CDKN1B
GO.0097305	Response to alcohol	5	0.0277	CCNA2, CCND1, CDK2, CDK4, CDKN1B
GO.1901991	Negative regulation of mitotic cell cycle phase transition	4	0.0287	CCNA2, CCND1, CDK2, CDKN1B
GO.0000082	G1/S transition of mitotic cell cycle	4	0.0326	CCND1, CDK2, CDK4, CDKN1B
GO.0048513	Organ development	11	0.0345	CCNA2, CCND1, CDK2, CDK4, CDKN1B, DLX3, KRT1, KRT10, PLCL2, PPP1R13L, SPINK5
GO.0060429	Epithelium development	7	0.0457	CCND1, CST6, DLX3, FABP5, KRT10, PTK6, SPINK5
GO.0048545	Response to steroid hormone	5	0.0482	CCNA2, CCND1, CDK2, CDK4, CDKN1B

[i] GO, Gene Ontology.

Performance of the biomarkers in clinical subtypes

Among the clinical factors, pathological T stage (T1 and T2-4), cancer status (with tumor and tumor free) and age (>57 years vs. ≤57 years) were all significantly associated with the recurrence status of patients with skin cancer (Fig. 6). Therefore, these 20 genes should to be considered in different clinical types to determine which clinical state they are more suitable for.

Figure 6 Survival curves of patients. Survival curves of the patients based on (A and B) cancer status, (C and D) pathological T-stage or (E and F) age, stratified by risk. HR, hazard ratio.

The results in Fig. 6 shows the predictive effect of these 20 gene biomarkers. Patients who were tumor free had improved survival compared with those with tumor. Similarly, patients with a T-stage of T1-3 exhibited improved survival compares with patients classed as T4, and patients ≤57 years old had improved survival compared with patients >57 years old.

Discussion

In the present study, a variety of statistical analysis methods were used (including LASSO regression, single-factor survival analysis, multi-factor Cox proportional hazards regression model and ROC curve analysis) to identify differences in the gene expression profiles of patients with skin cancer. A supervised cluster analysis method was used and 20 prognostic genes from the training samples were identified and verified against a test sample. The results showed that these 20 genes can stratify patients with skin cancer as high-risk and low-risk, which shows the feasibility of the mining method used in the present study.

From a biological point of view, the 20 prognostic genes identified in the present study successfully divided the patients with cancer according to the risk of recurrence, which may have important reference value for the treatment of recurrence of metastatic skin cancer, and may influence future clinical studies of skin cancer and drug development.

Among the 20 prognostic genes identified in the present study, several are closely associated with skin cancer and metastatic skin cancer, including DLX3, PTK6 and CST6 genes. For example, physical interaction between DLX3 and P53 on P21 promoter can enhance the expression of P21(53). Increasing DLX3 expression in keratinocytes produces a G1-S blockade associated with P53 signature transcriptional profiles (54).

Deletion of DLX3 promotes a mitogenic phenotype associated with constitutive activation of ERK (55). The loss of DLX3 expression in human skin cancer suggests that a DLX3-P53 interaction may be a primary regulatory axis of epidermis differentiation and that DLX3 may be a regulator of skin cancer development. Protein tyrosine kinase 6 (PTK6) is expressed in ~70% of cases of triple-negative breast cancer, in which it serves an important role in promoting metastatic lung colonization and survival (56). PTK6 inhibits the inhibition of metastasis of triple-negative breast cancer via SNAIL-dependent regulation of E-cadherin expression (57). Epigenetic changes associated with upregulation of CST6 gene expression may be accompanied by metastatic diffusion of primary tumor sites, and current studies have shown that methylation-dependent epigenetic silencing of CST6 represents an important mechanism for loss of CST6 during the development and/or progression to metastasis (58). Other gene prognostic biomarkers identified in the present study have potential research value and need further exploration and research.

Among these, several genes, including PLCL2, serine protease inhibitor Kazal type-5 (SPINK5) and KRT10 are associated with skin diseases. Inosine strongly enhances proliferation of human C32 melanoma cells through the PLCL2 and PI3K pathways (59). SPINK5 serves a crucial role in the timing of desquamation of the skin (60). Biallelic KRT10 mutations result in skin fragility caused by self-improving epidermolytic ichthyosis (61). Several other genes, such as transmembrane protein 45B (TMEM45B), CDKN1B, CKMT1A, FABP5 and PPP1R13L have also been shown to be associated with several types of cancer. TMEM45B has been shown to be abnormally expressed in gastric tumors and serves an important role in gastric tumorigenesis (62). The CCND1-A870G and CDKN1B-C79T polymorphisms are associated with breast cancer risk (63). The n335586/miR-924/CKMT1A axis contributes to migration and invasion of hepatocellular carcinoma cells (64). FABP5 serves an important role in the carcinogenesis and metastasis of cervical cancer, and may be a novel predictor for prognostic assessment of patients with cervical cancer (65). Polymorphisms of ERCC1, CD3EAP and PPP1R13L in chromosomal region 19q13.2-3 have previously been shown to exhibit a synergistic effect on apoptosis and DNA repair pathways (66).

In summary, the 20 gene biomarkers identified based on the LASSO algorithm can effectively predict the risk of patients with skin cancer and may be more convenient as a model of prognosis.

Acknowledgements

Not applicable.

Funding

This work was supported by the Key Project of China Ministry of Education for Philosophy and Social Science: Big Data Driven Risk Research on City's Public Safety (grant no. 16JZD023) and the Fundamental Research Funds for the Central Universities: Big Data Driven Risk Pre-Warning Research on City's Public Safety (grant no. 17LZUJBWZD012).

Availability of data and materials

The datasets used and/or analyzed during the present study are available from the corresponding author on reasonable request.

Authors' contributions

GL analyzed and interpreted the patient data was a major contributor in writing the manuscript. CL analyzed the data. HZ made substantial contributions to acquisition of data. ZZ was a major contributor in interpreting the data and in writing the manuscript. YS made substantial contributions to analysis and interpretation of data. All authors have read and approved the final manuscript.

Ethics approval and consent to participate

Not applicable.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References