Introduction
Cutaneous melanoma (CM) is an aggressive skin malignancy. In 2020, approximately 325,000 new melanoma cases and 57,000 melanoma-related deaths were estimated globally.1 A wide range of risk factors have been identified for CM, including cumulative exposure to ultraviolet (UV) radiation from sunlight, fair skin phenotype, atypical nevi, family history of CM, and immunosuppression. Given the ubiquity of UV exposure in daily life, CM can develop with relative ease in susceptible populations.2 Therefore, CM represents a substantial and considerable disease burden worldwide.
Melanoma cells exhibit extensive reprogramming of gene expression networks, driven by aberrant activation of oncogenic signaling pathways such as MAPK and PI3K/AKT.3 DNA methylation critically regulates gene expression. Concurrently, DNA methylation contributes to CM initiation and progression, as exemplified by promoter hypermethylation of the tumor suppressor gene PTEN.4
Psoriasis represents a chronic immune-mediated inflammatory skin disease.5 Similarly, psoriasis is also regulated by transcriptional and epigenetic mechanisms. Transcriptomic studies have identified upregulated genes in psoriatic lesions, including KRT16 and S100A7, while DNA methylation analyses have revealed epigenetic dysregulation of genes such as GJB2 and IRS1.6–9 Several studies have reported an increased overall cancer burden among patients with psoriasis;9,10 epidemiological evidence also suggests that psoriasis may be associated with a higher risk of CM,11,12 while the Mendelian randomization study by Li et al reached the opposite conclusion;13 whether the risk of CM development is associated with the presence of psoriasis remains controversial.
Molecular studies may shed light on the association between psoriasis and CM pathogenesis. As described above, the pathogenesis of both CM and psoriasis is closely associated with DNA methylation and gene expression. Previous studies also identified shared differentially expressed genes (DEGs) and shared differentially methylated genes (DMGs) between psoriasis and CM, which are implicated chemokine signaling and pro-inflammatory pathways.14,15 Similarly, experimental evidence indicates that the CXCL9–CXCR3 axis and JAK/STAT signaling may constitute a shared immune microenvironment between the two diseases.16 Beyond these, methylation-regulated genes (MRGs) which are characterized by coupled alterations in DNA methylation and gene expression, functioning as critical effectors that translate epigenetic dysregulation into disease-relevant phenotypes.17,18 Representative MRGs include MAGE family members in CM, and S100A9 in psoriasis, all of which are implicated in cell proliferation and immune-related pathways.4,19 Furthermore, Psoriatic lesions exhibit reproducible DNA methylation abnormalities,18 and resolved psoriatic lesions can retain disease-related transcriptional features, indicating the concept of molecular memory in the skin.20 Given that the MRGs, shared DEGs and DMGs implicated in both diseases are involved in common pathways including immune and inflammatory signaling, common MRGs (co-MRGs) may exist, whose DNA methylation-mediated dysregulation could link psoriasis to CM development via these pathways.4,14–20 However, most studies addressing psoriasis-related cancer risk have relied on single-omic expression analysis.21,22 Shared MRGs, defined by disease-associated transcriptional change together with regional DNA methylation alteration, remain insufficiently characterized in psoriasis and CM.
Moreover, compared with differential gene analysis alone, further exploration of MRGs may help provide biological plausibility for transcriptional alterations and offer more potential targets for intervention at the epigenetic level.14,23 Public multi-omics databases, including the Gene Expression Omnibus (GEO), ArrayExpress, The Cancer Genome Atlas (TCGA), and the Genotype-Tissue Expression (GTEx), provide readily available data resources for such investigations. We therefore hypothesized that co-MRGs exist between the two diseases, and to investigate this hypothesis based on the aforementioned databases. As two highly prevalent skin disorders with substantial health burdens, clarifying their molecular association can help identify high-risk populations and inform prevention and treatment.
Materials and Methods
Analysis Data Acquisition
Transcriptomic and epigenomic datasets were retrieved from the GEO (https://www.ncbi.nlm.nih.gov/geo/) and ArrayExpress (https://www.ebi.ac.uk/biostudies/arrayexpress) datasets. The search strategy employed the following terms: (Melanoma [MeSH Terms]) AND (Psoriasis [MeSH Terms]). This search was applied identically to both databases to identify relevant datasets.
RNA expression datasets included: GSE161683 (Psoriasis): a microarray dataset comprising 9 psoriatic and 9 healthy control skin samples; GSE100050 (CM): a microarray dataset comprising 6 CM and 6 normal skin samples. Methylation datasets included: E-GEOD-31835 (Psoriasis): a microarray-based methylation profiling dataset containing 19 psoriatic and 18 control samples; E-MTAB-14045 (CM): a microarray-based methylation profiling dataset containing 19 CM and 22 control samples; these datasets were incorporated into the main dataset. For external validation dataset of gene expression, RNA-sequencing (RNA-seq) data were obtained fromTCGA (https://portal.gdc.cancer.gov) and the GTEx (https://gtexportal.org) project and subsequently merged into a single dataset.
The TCGA melanoma cohort included 472 CM tumor samples and 1 normal sample, while GTEx provided 325 healthy skin control samples. Detailed information for each dataset is presented in Table S1. An overview of the data analysis pipeline is presented in Figure 1.Figure 1The workflow of current study.A workflow diagram linking methylation and transcriptome datasets to gene and pathway analyses.Two parallel workflow branches start with dataset sources and converge on downstream analyses.
Left branch, methylation dataset: E-MTAB-14045 and E-GEOD-31835 feed into, Differential analysis, then to, DMCs. Next is, Differentially methylated region (DMR) analysis, leading to, Differentially DMRs. This connects downward to, Methylation-regulated genes (MRGs), then to, Co-MRGs, then to, Over representation enrichment, Network analysis, Drug sensitivity analysis.
A side arrow from, Methylation-regulated genes (MRGs), points left to, Enrichment analysis (GO, KEGG). Right branch, transcriptome dataset: GSE100050 and GSE161683 feed into, Differential analysis, then to, Overlapped genes. A connector from this branch leads to, Differentially expressed genes (DEGs), which then connects to, Co-MRGs.
From, Co-MRGs, a vertical sequence of analysis boxes is shown to the right: GSEA analysis; Expression and identification analysis; Immune infiltration analysis; Survival analysis; Immunohistochemistry; Cell line expression profiling. A bracket groups the first three items, GSEA analysis, Expression and identification analysis and Immune infiltration analysis and a right-side arrow points into this grouped section from two stacked boxes labeled, gtexKallistoestcounts and, TCGA-SKCM.starcounts. A separate arrow from the same right-side source points to, Survival analysis.
Data Preprocessing
The microarray data were pre-filtered before extraction. All RNA microarray datasets were pre-normalized by the original data submitters (Log2-RMA). Methylation M-values were either directly obtained from the datasets or calculated from beta values (Figure S1).
TCGA and GTEx RNA-seq data were processed using the edgeR R package to filter low-expression genes and for normalization. Low-expression genes were filtered using the filterByExpr function with default settings. The filtered count matrix was then normalized using the cpm() function, which converts counts to log2-CPM values based on effective library sizes.24 To ensure comparability, the merged RNA-seq datasets were corrected for batch effects using the ComBat algorithm from the sva R package.25 Data source (TCGA vs GTEx) was specified as the batch variable.
The effectiveness of normalization and batch correction is demonstrated in Figure S2.
Identification of Differentially Expressed Genes (DEGs) and Differentially Methylated Positions (DMPs)
Principal Component Analysis (PCA) was initially performed to visualize the overall variance in gene expression and methylation patterns between case and control groups. DEGs and DMPs were identified using the limma R package.26 The significance thresholds were set as an adjusted p-value (False Discovery Rate, FDR) < 0.05 and |log2 fold change (FC)| ≥ 1.
Gene Function Enrichment Analysis
The Gene Ontology (GO) functional enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were performed on the MRGs, while Gene Set Enrichment Analysis (GSEA) was conducted based on the logFC of all genes, using the clusterProfiler R package.27 GO and KEGG analyses were performed separately on the DEGs identified from the RNA microarray data of the psoriasis dataset and the melanoma dataset. The threshold for significance was set at q-value < 0.05. GSEA was additionally conducted on the merged TCGA-GTEx dataset to characterize the up- and down-regulation patterns of pathways.
Significance was set at FDR < 0.05.
Differentially Methylated Regions (DMRs) Annotation and Identification of MRGs
To further refine the methylation analysis, the DMRcate R package was used to identify DMRs between case and control groups, with a FDR threshold of 0.05.28 DMR calling thresholds were adapted to the CpG density of each platform. For the E-GEOD-31835 (27K array), in which CpG probes are relatively sparsely distributed, DMRs were defined as regions containing at least two significant DMPs within a 1000 bp window to minimize the risk of missing potential differential DMRs. For the E-MTAB-14045 dataset (450K array), which features a substantially higher density of CpG sites and thus a greater likelihood of observing significant DMPs by chance, a more stringent threshold of at least five significant DMPs within the same window size was applied to maintain comparable specificity and avoid false-positive inflation.
These DMRs were then annotated to genes using the hg19/GRCh37 database. Genes that were both DEGs and overlapped with annotated DMRs were identified as MRGs for each disease. The intersection of MRGs between psoriasis and CM was defined as co-MRGs.
Network Analysis
Network analyses were conducted based on the identified co-MRGs. The correlation of co-MRGs expression was analyzed using the R package “circlize” and visualized via a chord diagram.29 A protein-protein interaction (PPI) network was constructed using the STRING database (https://string-db.org/) with a minimum required interaction score of 0.4. The resulting network was imported into Cytoscape for visualization and analysis.30 To delineate the regulatory relationships between the co-MRGs and miRNAs, potential interactions were retrieved from TarBase databases to construct the miRNA-mRNA regulatory network.31 In parallel, a functional association network for the co-MRGs was generated using the GeneMANIA database (https://genemania.org/).
For drug target analysis, the identified co-MRGs were submitted to the Drug-Gene Interaction Database (DGIdb; https://www.dgidb.org/) to retrieve information on existing drugs or compounds that target these genes. Subsequently, we conducted an additional drug-sensitivity analysis of these co-MRGs in melanoma using the Comprehensive Pancancer Analysis of Drug Sensitivity (CPADS) web tool, stratified by high and low expression groups, with data sourced from the Genomics of Drug Sensitivity in Cancer (GDSC). The IC50 value was used as a measure of drug resistance.
Expression and Diagnostic Performance Analysis
The expression of co-MRGs was further validated in the external TCGA-GTEx dataset. Their diagnostic performance was assessed by plotting Receiver Operating Characteristic (ROC) curves using the pROC R package and calculating the Area Under the ROC Curve (AUC), where a value closer to 1 indicates superior diagnostic ability.32 Bar plots were generated to illustrate the expression levels of these genes between case and control.
Immune Cell Infiltration Analysis
Using the merged TCGA-GTEx dataset, we quantified the relative abundances of 22 tumor-infiltrating immune cell types by the CIBERSORT algorithm, as implemented in the CIBERSORT R package, and the ssGSEA method, as implemented in the GSVA R package.33,34 Furthermore, the correlations between the expression levels of the identified co-MRGs and the degree of immune cell infiltration were analyzed.
Survival Analysis
For survival analysis, 439 CM patients with available survival data from the TCGA cohort were stratified into high- and low-expression groups based on the median expression level of each validated co-MRG. Survival curves were generated using the Kaplan-Meier method with the survival R package.35 Differences in survival distributions between groups were compared using the Log rank test.
Immunohistochemistry (IHC) and Cell Line Expression Profiling
To validate the protein expression patterns of our candidate genes at both the tissue and cellular levels, we utilized the Human Protein Atlas (HPA, https://www.proteinatlas.org), a comprehensive open-access resource integrating antibody-based IHC and transcriptomics data. HPA’s cell line clustering analysis was performed, which groups genes using the Louvain algorithm based on their expression similarity across the full panel of cancer cell lines profiled by HPA. The resulting clusters were visualized via Uniform Manifold Approximation and Projection (UMAP) dimensionality reduction, with GO functional annotation to characterize the biological processes associated with each cluster.
This approach allowed us to contextualize the expression patterns of our candidate genes within a malignant cellular context.
Statistical Analysis
All statistical analyses were performed using R software (version 4.4.2). Comparisons between two groups were assessed using the Wilcoxon rank-sum test or Student’s t-test, as appropriate based on data distribution. Correlation analyses were performed using Spearman’s rank correlation coefficient.
A two-sided p-value of less than 0.05 was considered statistically significant.
Results
Identification of DEGs and DMPs
PCA initially revealed differences in expression and methylation levels between case and control groups across the four datasets (Figure S3). In the RNA datasets, a total of 598 DEGs were identified in the psoriasis dataset GSE161683, including 340 upregulated and 258 downregulated genes (Figure 2a and b). In the CM dataset GSE100050, 1821 DEGs were identified, of which 1146 were upregulated and 675 were downregulated (Figure 2c and d).
Regarding the methylation datasets, the psoriasis dataset E-GEOD-31835 yielded 252 DMPs, comprising 90 hypermethylated and 162 hypomethylated sites (Figure 2e and f). The CM dataset E-MTAB-14045 identified 443 DMPs, including 360 hypermethylated and 83 hypomethylated sites (Figure 2g and h).Figure 2Differential expression and differential methylation analyses in main dataset. Volcano plots and heatmaps of differentially expressed genes (DEGs) in (a and b) psoriasis and (c and d) cutaneous melanoma (CM), respectively.
Volcano plots and heatmaps of differentially methylated positions in (e and f) psoriasis and (g and h) CM, respectively. (i) Venn diagram showing the overlap between DEGs and differentially methylated regions (DMRs). DMRs in the plot represent genes annotated to differentially methylated regions.Different types of data visualizations with four volcano plots, four heatmaps and one Venn diagram.A) Volcano plot titled, GSE161683 (Psoriasis, RNA dataset).
X axis label: log2FoldChange, unit not shown; range approximately negative 3 to 3. Y axis label: minus log10Pvalue, unit not shown; range approximately 0 to 15. Legend title: State, entries Down, Stable, Up.
Two vertical dashed thresholds bracket x equals 0 and one horizontal dashed threshold is near y equals 5. Points are densest near x equals 0 at low y, with many points beyond the thresholds on both negative and positive x reaching up to about y equals 15. B) Heatmap under the same dataset heading as A.
X axis label not shown; unit not shown. Y axis label not shown; unit not shown. Left-side dendrogram.
Top annotation labeled, Group, with CASE and NG. Right-side heatmap scale shows values from negative 2 to 2, unit not shown. Columns split into two main groups aligned with the Group bar, with large row blocks showing opposite high versus low patterns between the two column groups.
C) Volcano plot titled, GSE100050 (CM, RNA dataset). Same axis labels and units as A. X axis range approximately negative 5 to 5.
Y axis range approximately 0 to 15. Same State legend and dashed thresholds as A. Dense center near x equals 0 with broad spread; many points on both sides exceed the horizontal threshold, with the highest points near y equals 15.
D) Heatmap under the same dataset heading as C. Same layout as B, including Group with CASE and NG and a scale from negative 2 to 2, unit not shown. Compared with B, the dominant block structure spans more rows, with a large region where one column group is consistently high while the other is consistently low.
E) Volcano plot titled, E-GEOD-31835 (Psoriasis, Methylation dataset). X axis label: log2FoldChange, unit not shown; range approximately negative 3 to 3. Y axis label: minus log10Pvalue, unit not shown; range approximately 0 to 6.
Same State legend as A. Two vertical dashed thresholds bracket x equals 0 and one horizontal dashed threshold is near y equals 1. Most points cluster near x equals 0 at low y, with smaller sets on both sides rising to about y equals 6.
F) Heatmap under the same dataset heading as E. Same layout elements as B, including Group with CASE and NG and a scale from negative 2 to 2, unit not shown. Compared with B, the heatmap shows more alternating bands across rows and a less uniform block separation between the two column groups.
G) Volcano plot titled, E-MTAB-14045 (CM, Methylation dataset). Same axis labels and units as E. X axis range approximately negative 4 to 4.
Y axis range approximately 0 to 15. Same State legend as A. Two vertical dashed thresholds bracket x equals 0 and one horizontal dashed threshold is near y equals 1.
Point cloud forms a pronounced U-shaped envelope, with lowest y near x equals 0 and higher y toward both extremes; the positive x side has a dense set reaching near y equals 15. H) Heatmap under the same dataset heading as G. Same layout elements as B, including Group with CASE and NG and a scale from negative 2 to 2, unit not shown.
Compared with F, separation between the two column groups is stronger, with a large row block showing consistently opposite values between the groups. I) Venn diagram with four sets labeled DEGs-Psoriasis, DEGs-CM, DMRs-CM and DMRs-Psoriasis. Counts shown: DEGs-Psoriasis only 174; DEGs-CM only 891; DMRs-CM only 11879; DMRs-Psoriasis only 56.
Pair overlaps shown: DEGs-Psoriasis with DEGs-CM 88; DEGs-CM with DMRs-CM 705; DMRs-CM with DMRs-Psoriasis 182; DEGs-Psoriasis with DMRs-Psoriasis 7. Additional overlap counts shown inside intersections are 93, 223, 25, 13, 3, 3 and the central overlap is 3.
Identification of DMRs and Co-MRGs
Using DMRcate, DMRs associated with psoriasis and CM were further identified from the DNA methylation microarray dataset in the main dataset. A total of 251 DMRs associated with psoriasis and 16,569 DMRs associated with CM were detected, which were annotated to 296 and 13117 genes, respectively (Tables S2–S3). The intersections of DEGs and DMRs, designated as MRGs, comprised 20 and 826 genes for psoriasis and CM, respectively, and are presented in a Venn diagram (Figure 2i) and Tables S4-S5.
Ultimately, three co-MRGs (KRT6B, TSPAN8, CHI3L2) shared between psoriasis and CM were identified (Figure 2i and Table 1).Table 1Identification of Co-MRGs in Psoriasis and Cutaneous Melanoma from the Main DatasetGeneLogFCAveExpr (log2-RMA)tAdj.p.valUp- or Down RegulationHyper-or Hypo MethylationPsoriasisCHI3L22.428.826.540.00UpHypoKRT6B2.3012.099.110.00UpHypoTSPAN8−1.825.93−6.720.00DownHyperCutaneous melanomaTSPAN8−2.809.82−6.680.00DownHyperCHI3L22.547.344.810.01UpHypoKRT6B−2.0710.90−4.560.01DownHyper
Enrichment Analysis
GO and KEGG functional analyses were performed on the MRGs of psoriasis and CM, respectively. For psoriasis, the enrichment analysis mainly identified immune- and inflammation-related terms, such as RAGE receptor binding and Toll-like receptor binding, as well as terms associated with epidermal dysfunction, including structural constituent of skin epidermis and cornified envelope formation (Figure S4a). For CM, the GO enrichment analysis primarily yielded terms related to the extracellular matrix, and also identified immune regulation and skin development-related terms, such as leukocyte cell–cell adhesion and epidermis development (Figure S4b).
KEGG pathway analysis revealed pathways highly associated with tumor occurrence and metastasis, including central carbon metabolism in cancer and cell adhesion molecule (CAM) interactions (Figure S4c). In addition, an independent GO enrichment analysis was performed specifically on the three co-MRGs, and the top 10 most significant terms were visualized using a mulberry plot to illustrate the core functions of these genes (Figure S4d). KRT6B was primarily involved in skin structure and differentiation, CHI3L2 participated in glucose metabolism, and TSPAN8 was associated with hemostasis regulation and wound healing.
Network Analysis
We constructed expression correlation networks for the co-MRGs in the two RNA microarray datasets (Figure 3a and b). The results showed inconsistent expression patterns across the datasets. In psoriasis, KRT6B expression was positively correlated with CHI3L2 and negatively correlated with TSPAN8.
In contrast, the opposite correlations were observed in CM. TSPAN8 and CHI3L2 showed a negative correlation in both psoriasis and CM. The PPI network revealed that the three co-MRGs belonged to distinct protein interaction networks (Figure 3c); however, functional association analysis identified a certain degree of co-expression network among them (Figure 3d).
We further constructed miRNA–mRNA regulatory networks for the co-MRGs and found that only hsa-let-7i-5p and hsa-let-7g-5p were capable of regulating both KRT6B and CHI3L2. (Figure 3e). Drug target analysis identified two FDA-approved drugs targeting CHI3L2 (dipyridamole and ticagrelor), both associated with antithrombotic effects,36 as well as two equilibrative nucleoside transporter (ENT) inhibitors with potential anti-inflammatory effects, NBTGR and KF2434537,38 (Figure 3f).
The three drugs with the most significant differences in drug resistance between the high- and low-expression subgroups of each co-MRG are presented in Figure S5. Notably, methotrexate, a drug used for the treatment of severe and recalcitrant psoriasis,39 demonstrated the most significant differential response between high- and low-CHI3L2 expression groups, with higher CHI3L2 expression associated with markedly decreased drug sensitivity (p<0.001).Figure 3Network analyses. (a and b) Correlation analysis of common methylation-regulated genes (co-MRGs) expression in psoriasis and cutaneous melanoma, respectively.
(c) Protein-protein interaction network of co-MRGs. (d) Functional association network of co-MRGs. (e) miRNA-mRNA regulatory network of co-MRGs.
(f) Drug-target network.Six-part network analysis: two chord correlation plots and four interaction networks for co-MRGs.A) Circular chord correlation plot for three labeled genes, KRT6B, CHI3L2 and TSPAN8. A vertical color key labeled, Correlation, shows a scale from 1 at the top through 0 to minus 1 at the bottom. Three wide chords connect the three gene sectors, with one chord shaded toward the positive end of the key and two chords shaded toward the negative end of the key.
B) Circular chord correlation plot for the same three labeled genes, KRT6B, CHI3L2 and TSPAN8, with the same vertical color key labeled, Correlation, scaled from 1 through 0 to minus 1. Three wide chords connect the three gene sectors, with two chords shaded toward the positive end of the key and one chord shaded toward the negative end of the key. C) Protein protein interaction style network with many circular nodes labeled with gene names.
A dense central cluster contains multiple SPRR and KRT labeled nodes connected by many lines. Three co-MRG nodes are shown as labeled circles: KRT6B and KRT6C are positioned at the edge of the dense cluster, while TSPAN8 is separated to the right and connects by a single line to ZBED3. CHI3L2 is separated toward the lower left and connects by a single line to DENND2D.
D) Functional association network with many oval nodes connected by thin lines. Three larger labeled nodes, KRT6B, CHI3L2 and TSPAN8, sit within a broader web of smaller unlabeled or less legible nodes. KRT6B is positioned left of center with multiple connections into the network.
CHI3L2 is positioned right of center with multiple connections into a right side cluster. TSPAN8 is positioned upper right with fewer connections, including a long connection extending to a small node at the far right. E) miRNA mRNA regulatory network drawn as a layered graph.
Two larger nodes at the top connect downward through many directed links to numerous small circular nodes arranged in rows. Many of the small nodes have text labels that are not legible at this size. The overall structure shows extensive fan out from the top nodes to many downstream nodes.
F) Drug target network centered on a node labeled, CHI3L2, with multiple spokes to surrounding nodes. Surrounding nodes include labels, dipyridamole, ticagrelor, NBTGR and KF24345, plus additional drug name labels that are not legible. A legend titled, Node types, shows three marker categories: Gene, Approved drug and Not approved drug.
Validation Analysis
A combined dataset from TCGA and GTEx was used as the validation cohort, comprising 472 CM cases and 326 normal controls. The results showed that KRT6B and TSPAN8 remained significantly downregulated (FDR<0.001), and CHI3L2 remained significantly upregulated (FDR<0.001) (Table 2). GSEA was performed on the validation dataset.
Pathways involving KRT6B were significantly upregulated and associated with epidermal and skin development as well as keratinization, and pathways involving CHI3L2, including carbohydrate metabolism, were also upregulated; in contrast, the integrin-binding pathway involving TSPAN8 was significantly downregulated (Figure 4a-c). Diagnostic analysis revealed that KRT6B (AUC = 0.822) and CHI3L2 (AUC = 0.867) exhibited superior diagnostic performance for CM compared to TSPAN8 (AUC = 0.606) (Figure 4d-f). The significantly reduced expression of KRT6B and TSPAN8, and the significantly elevated expression of CHI3L2 in cases are shown in Figure 4g-i.Table 2Expression Levels of Co-MRGs in the Validation Dataset (TCGA-GTEx)logFCAveExpr (log2- RMA)tadj.P.ValCHI3L22.774.80252.070.00KRT6B−2.336.05−66.570.00TSPAN8−1.205.00−20.490.00Figure 4TCGA–GTEx validation analysis.
GSEA plots for (a) KRT6B, (b) TSPAN8, and (c) CHI3L2. Diagnostic performance of (d) KRT6B, (e) TSPAN8, and (f) CHI3L2 in cutaneous melanoma (CM). Expression of (g) KRT6B, (h) TSPAN8, and (i) CHI3L2 in the validation cohort.Multiple plots of GSEA curves, ROC curves and gene expression box plots for KRT6B, TSPAN8 and CHI3L2.A) GSEA line graph with three stacked sections.
X-axis label: Rank in Ordered Dataset, unit not shown. Y-axis label: Running Enrichment Score, unit not shown, with labeled values 0.25, 0.50, 0.75. Four enrichment curves start near 0.75 to 0.85 at the left edge and decrease toward about 0.00 at the far right.
Pathway labels: epidermal cell differentiation, epidermis development, keratinization, keratinocyte differentiation, skin development. Inset table columns: pvalue and p.adjust. Five rows show pvalue equals 1e minus 10 and p.adjust equals 1.454e minus 08.
Middle section shows multiple rows of vertical tick marks across the rank axis. Bottom section label: Ranked List Metric, with a distribution high at the left and tapering toward the right. B) GSEA line graph with the same three stacked sections and the same x-axis label as A).
Y-axis label: Running Enrichment Score, unit not shown, with labeled values negative 0.4, negative 0.2, 0.0, 0.2. One enrichment curve starts near 0.2, drops below 0.0 to about negative 0.4 in the mid to right portion, then rises near the far right to above 0.0. Pathway label: integrin binding.
Inset table shows pvalue equals 0.004936 and p.adjust equals 0.04655. Middle section shows a dense band of vertical tick marks. Bottom section label: Ranked List Metric.
C) GSEA line graph with the same three stacked sections and the same x-axis label as A). Y-axis label: Running Enrichment Score, unit not shown, with labeled values 0.0, 0.2, 0.4, 0.6. Two enrichment curves rise from near 0.0 to about 0.6 near the left to mid portion, then decline toward about 0.1 by the far right.
Pathway labels: carbohydrate binding, carbohydrate derivative catabolic process. Inset table columns: pvalue and p.adjust. Row 1: pvalue equals 0.015442 and p.adjust equals 0.03534.
Row 2: pvalue equals 0.004067 and p.adjust equals 0.03534. Middle section shows two bands of vertical tick marks. Bottom section label: Ranked List Metric.
D) ROC curve line graph. X-axis label: Specificity, unit not shown, ranging from 1.0 to 0.0. Y-axis label: Sensitivity, unit not shown, ranging from 0.0 to 1.0.
A diagonal reference line is present. The ROC curve rises steeply near specificity close to 1.0 and then approaches the top boundary. Text: AUC equals 0.822.
E) ROC curve line graph with the same axes and ranges as D). The ROC curve increases gradually and stays closer to the diagonal reference line for much of the range. Text: AUC equals 0.606.
F) ROC curve line graph with the same axes and ranges as D). The ROC curve rises almost vertically near specificity close to 1.0 up to sensitivity near 0.9, then runs nearly flat near the top before reaching the upper right corner. Text: AUC equals 0.867.
G) Box plot comparing two groups. X-axis label: Group, unit not shown, categories control and case. Y-axis label: Log2CPM, unit not shown, labeled values 4, 5, 6, 7, 8.
Text above: p equals 6.3e minus 07. Control box: about 6.2 to 6.4 with median near 6.3 and whiskers about 5.8 to 6.6. Case box: about 5.0 to 6.7 with median near 6.0 and whiskers about 4.0 to 8.0.
H) Box plot comparing two groups with the same x-axis label and categories as G). Y-axis label: Log2CPM, unit not shown, labeled values 3, 4, 5, 6, 7. Text above: p equals 0.00024.
Control box: about 5.0 to 5.6 with median near 5.4 and whiskers about 4.6 to 6.4. Case box: about 4.3 to 5.6 with median near 5.0 and whiskers about 3.5 to 7.2. I) Box plot comparing two groups with the same x-axis label and categories as G).
Y-axis label: Log2CPM, unit not shown, labeled values 6, 7, 8, 9, 10. Text above: p less than 2e minus 16. Control box: about 6.3 to 7.0 with median near 6.6 and whiskers about 6.0 to 8.4.
Case box: about 7.1 to 9.1 with median near 8.2 and whiskers about 6.2 to 10.2.
Immune Infiltration Analysis
Immune infiltration in the validation datasets was assessed by ssGSEA and CIBERSORT. According to the ssGSEA infiltration scores, most immune cell types were significantly higher in CM cases (p<0.01), except for CD8⁺ T cells (Figure S6a). Moreover, these scores were mostly positively correlated (Figure S6b).
KRT6B expression was broadly negatively correlated with these scores, whereas TSPAN8 and CHI3L2 were positively correlated (p<0.05) (Figure S6c-e). To further validate these findings at the level of estimated cellular proportions, we applied the CIBERSORT algorithm (Figure 5) and confirmed that the proportions of M0 macrophages, M1 macrophages, plasma cells, and resting CD4⁺ memory T cells were significantly higher in CM cases (p<0.01) (Figure 5b). Among these, correlation analysis revealed that KRT6B was negatively correlated with the levels of M1 macrophages and resting CD4⁺ memory T cells (p<0.01), while TSPAN8 was positively correlated with these two cell types (p<0.01) but negatively correlated with the levels of M0 macrophages and plasma cells (p<0.001).
CHI3L2 was positively correlated with the levels of M0 macrophages and resting CD4⁺ memory T cells (p<0.05) (Figure 5d–f).Figure 5CIBERSORT-based immune infiltration analysis in the CM-related validation dataset. (a) Proportions of infiltrating immune cells in each sample. (b) Comparison of immune cell composition between the case and control groups.
(c) Correlation matrix of immune cell infiltration proportions. Correlation between immune cell infiltration and the expression of (d) KRT6B, (e) TSPAN8, and (f) CHI3L2. * p <0.05, ** p <0.01, *** p <0.001.Different plots show immune cell composition, group comparisons, correlations and gene associations.A) Stacked bar chart titled, Immune Cell Composition. Two grouped sections labeled, Case and, Normal.
The y axis label is, Proportion, with values from 0.00 to 1.00. The x axis label is not legible, showing many narrow stacked bars representing samples. A legend titled, Cell Type, lists: B cells memory, B cells naive, Dendritic cells activated, Dendritic cells resting, Eosinophils, Macrophages M0, Macrophages M1, Macrophages M2, Mast cells activated, Mast cells resting, Monocytes, Neutrophils, Natural killer cells activated, Natural killer cells resting, Plasma cells, T cells CD4 memory activated, T cells CD4 memory resting, T cells CD4 naive, T cells follicular helper, T cells gamma delta, T cells regulatory Tregs, T cells CD8 and T cells receptor gamma delta.
Each sample bar is subdivided into segments whose combined height equals 1.00. B) Grouped dot and bar plot comparing immune cell fractions by group. Legend label, group, with categories, Case and, Normal.
The x axis label is, Cell Type, with the same immune cell categories as A). The y axis label is, Fraction, ranging from 0.0 to 0.4. For each cell type, two vertical point distributions are shown with a central bar.
Significance markers appear above several cell types as single, double and triple asterisks. C) Correlation heatmap matrix of immune cell infiltration proportions. Both axes list the same immune cell types as in A).
A vertical scale at the right is labeled from negative 1 through 0 to 1. Each cell contains a numeric correlation value that is not legible. The diagonal shows self correlations from top left to bottom right.
D) Lollipop style correlation plot titled, Spearman Correlation with KRT6B. Legend titled, Correlation, with categories, Negative and, Positive. The x axis label is, Correlation Coefficient, ranging from negative 0.50 to 0.50.
The y axis label is, Variables, listing immune cell types including: Dendritic cells activated, Mast cells resting, Mast cells activated, Neutrophils, Dendritic cells resting, T cells CD4 naive, T cells gamma delta, B cells memory, T cells CD4 memory activated, Natural killer cells resting, T cells regulatory Tregs, Monocytes, Macrophages M0, T cells follicular helper, Natural killer cells activated, Plasma cells, B cells naive, T cells CD8, T cells CD4 memory resting, Eosinophils, Macrophages M1 and Macrophages M2. Each variable has a point connected to the zero line by a horizontal segment, positioned left for negative and right for positive, with asterisks shown for some variables. E) Lollipop style correlation plot titled, Spearman Correlation with TSPAN8.
Same layout, legend categories, x axis label and x axis range as D). The y axis label is, Variables, listing immune cell types including: Dendritic cells resting, Macrophages M1, Mast cells resting, T cells CD4 memory resting, B cells naive, Natural killer cells resting, B cells memory, T cells CD4 memory activated, T cells CD8, Macrophages M2, Neutrophils, T cells gamma delta, T cells follicular helper, T cells CD4 naive, Eosinophils, Mast cells activated, Natural killer cells activated, Dendritic cells activated, T cells regulatory Tregs, Plasma cells, Monocytes and Macrophages M0. Points and horizontal segments indicate negative or positive correlations relative to zero, with asterisks shown for some variables.
F) Lollipop style correlation plot titled, Spearman Correlation with CHI3L2. Same layout, legend categories, x axis label and x axis range as D). The y axis label is, Variables, listing immune cell types including: Macrophages M0, T cells CD4 memory resting, B cells memory, Mast cells resting, Mast cells resting, B cells naive, Macrophages M1, Macrophages M2, Monocytes, Dendritic cells resting, Neutrophils, Natural killer cells activated, T cells gamma delta, Plasma cells, T cells regulatory Tregs, T cells follicular helper, T cells CD4 naive, Dendritic cells activated, Eosinophils, Mast cells activated, T cells CD8, Natural killer cells resting and T cells CD4 memory activated.
Each variable is plotted as a point with a horizontal segment to the zero line, with asterisks shown for some variables.
Survival Analysis
Based on the follow-up data of 439 cases from TCGA, survival analysis was performed between the high- and low-expression groups of co-MRGs (Figure S7). Although none of the results reached statistical significance, high KRT6B expression showed a trend toward an association with worse survival in CM (log-rank p = 0.072), whereas the opposite results was observed for CHI3L2, with high expression showing a trend toward an association with better survival in CM (log-rank p = 0.081).
IHC Validation and Cell Line Expression Profiling Analysis based on HPA
KRT6B exhibited High staining in the corneal, granular, and spinous layer cells of normal skin, but was Not detected in basal layer cells or melanocytes (Figure 6a and b). We selected NOS (Not Otherwise Specified) CM as a surrogate validation tissue. KRT6B staining in CM tissues ranged from Not detected to Medium (Figure 6c and d).
In contrast, TSPAN8 was Not detected in both normal skin (Figure 6e and f) and CM tissues (Figure 6g and h). CHI3L2 was Not detected in normal skin (Figure 6i and j) but showed Medium to High staining in CM tissues (Figure 6k and l). Based on cell line clustering analysis from the HPA database, KRT6B was assigned to a gene cluster functionally annotated as “Squamous epithelial cell – Keratinization”; TSPAN8 was assigned to a cluster annotated as “Stomach and colon cancer – digestion”, and CHI3L2 to a cluster annotated as “T-cell receptor” (Figure 6m-o).Figure 6Immunohistochemistry (IHC) validation and cell line expression profiling analysis.
IHC staining of KRT6B (antibody HPA045697) in (a and b) normal skin and (c and d) cutaneous melanoma (CM) tissues. IHC staining of TSPAN8 (antibody HPA044337) in (e and f) normal skin and (g and h) CM tissues. IHC staining of CHI3L2 (antibody HPA005443) in (i and j) normal skin and (k and l) CM tissues.
Expression profiling of (m) KRT6B, (n) TSPAN8, and (o) CHI3L2 across different cell lines.Composite figure: 12 IHC tissue micrographs (a to l) plus 3 cell-line expression scatter plots (m to o).
Discussion
Integrated multi-omics analysis, by simultaneously interrogating multiple molecular layers, enables systematic dissection of cancer-driving regulatory mechanisms and enhances the sensitivity and specificity of biomarker discovery beyond single-omics approaches.40,41 CM represents a substantial global health burden among skin malignancies. However, the molecular mechanisms underlying its potential association with psoriasis remain poorly understood from a multi-omics perspective. In this study, we integrated transcriptomic and DNA methylation data to systematically investigate the shared molecular features between psoriasis and CM, and identified three co-MRGs: KRT6B, CHI3L2, and TSPAN8.
Notably, KRT6B and CHI3L2 exhibited robust diagnostic performance in CM, and their differential expression was further validated by IHC.
Enrichment results indicated epidermal differentiation and epithelial–immune communication as shared themes across diseases. Psoriasis MRGs were linked to inflammation and epidermal barrier, while CM MRGs were associated with matrix remodeling, adhesion, metabolism, and epidermal development. This pattern parallels the psoriasis keratinocyte–immune circuit and the CM microenvironment-remodeling models.42,43 Earlier reports also support convergent signatures in cell-cycle regulation and inflammatory signaling, consistent with our results.22,44
In psoriasis, KRT6B is upregulated along with KRT16 and KRT17 as a barrier alarmin, promoting keratinocyte hyperproliferation and Th1/Th17-driven inflammation.45 Contrary to our initial hypothesis, the opposing methylation and expression patterns of KRT6B between psoriasis and CM likely reflect the protective role of psoriasis against CM development, as some previous studies suggested.13 To date, no direct evidence has linked DNA methylation alterations to KRT6B expression in psoriasis. However, evidence from atopic dermatitis has revealed a comparable hypomethylation and concomitant upregulation of KRT6A, a highly homologous member of the KRT family;46 additionally, differential methylation of the KRT6B promoter has been documented in common warts.47 Similarly, studies on KRT6B in CM remain scarce, with only transcriptome-based database analyses reporting downregulation of KRT6B in CM, a finding consistent with our results.48 Although immune infiltration correlates of KRT6B in CM have not been investigated, a related keratin family member, KRT16, was found to be significantly negatively correlated with infiltration of CD4⁺ T cells and macrophages in CM.48 In addition, higher CD4⁺ T cell infiltration in CM was associated with better prognosis.49 We further observed a trend toward an association between high KRT6B expression and poor prognosis in CM, corroborating previous reports. Moreover, KRT6B expression was significantly higher in metastatic CM compared with non-metastatic CM.50,51 These findings implying that KRT6B may share similar immune-related patterns, whereby high KRT6B expression in psoriasis may suppress CD4⁺ T-cell infiltration, thereby contributing to poor prognosis in CM, although high KRT6B expression itself may be associated with a reduced risk of CM development.
The hypotheses require further validation in future studies.
CHI3L2 encodes a chitinase-like lectin that binds chitooligosaccharides with micromolar affinity.52 In present study, CHI3L2 represents a biologically plausible candidate, with hypomethylation and upregulation in both psoriasis and CM and strong diagnostic value. Although CHI3L2 is poorly characterized in both diseases, one study identified it as a core differential gene in psoriasis.53CHI3L2 encodes a chitinase-like protein implicated in tumor and immune contexts in several cancers. In invasive ductal breast carcinoma, CHI3L2 modulates tumor cell behavior via ERK1/2 and STAT3 signaling;54 in glioma, it serves as a prognostic biomarker linked to immune infiltration,55 supporting its broad relevance to tumor biology and immune regulation, though the functional significance of CHI3L2 in CM has yet to be determined.
We further identified positive associations between CHI3L2 and the infiltration of CD4⁺ T cells. In addition, while CHI3L2 is upregulated in glioma with adverse prognosis,55 its upregulation in CM showed a trend toward association with better prognosis in our study. Although no similar reports have been documented in the context of psoriasis or CM, a comparable pattern has been observed in studies on non-small cell lung cancer, where CD4⁺ T cells were associated with a better prognosis.56 The high-expression group showed diminished sensitivity to methotrexate, a conventional agent for severe psoriasis that has also been implicated as a potential risk factor for melanoma,57 suggesting that patients with severe psoriasis may need to be aware of the risk of developing CM when using this drug.
TSPAN8 is a membrane tetraspanin that organizes tetraspanin-enriched microdomains and modulates integrin signaling, adhesion, motility, invasion and extracellular-vesicle biology.58 No previous studies have reported an association between TSPAN8 and psoriasis; only one study identified a TSPAN8 mutation in keloid.59 Functional CM studies have reported that TSPAN8 promotes invasive behavior through β1-integrin-associated adhesion signaling, proMMP-9 activation and basement-membrane proteolysis.60,61 In human CM, TSPAN8 expression has also been associated with increased metastatic risk and cancer-related death.62 However, in our study, TSPAN8 was hypermethylated and downregulated in both disease contexts and remained downregulated in the TCGA–GTEx validation cohort. One possible explanation is that TSPAN8 expression varies across CM progression. It may be reduced in bulk tumor-versus-normal comparisons but become functionally relevant in invasive or metastatic subpopulations.
This discrepancy underscores the need to elucidate the context-dependent functions of TSPAN8 in the setting of psoriasis across different stages of tumor development.
Our findings further support a potential methylation–expression regulatory relationship between psoriasis and CM. We propose that the differential contribution of molecular alterations across distinct clinical contexts may determine the nature of the relationship between the two diseases. Specifically, when KRT6B hypomethylation and upregulation predominate, psoriasis may exert a protective effect against CM at the molecular level.
This dual risk–protective paradigm is supported by GWAS-based evidence on HLA-DOB, which confers psoriasis risk but serves as a protective factor against CM development.13 Conversely, when key molecules such as TSPAN8 and CHI3L2, which exhibit concordant epigenetic and transcriptional alterations in both psoriasis and CM, predominate and their effects outweigh those of KRT6B, psoriasis may instead serve as a risk factor for CM. Of note, the hypermethylation and concomitant downregulation of KRT6B and CHI3L2 have been previously documented in independent studies.63,64 Network analysis revealed that each of the three genes may be subject to its own relatively independent regulatory module. This proposition may partly account for the inconsistent conclusions regarding the comorbidity between psoriasis and CM in previous studies.
Our hypothesis is illustrated in Figure 7.Figure 7Schematic diagram of the proposed mechanism. We propose a mechanistic hypothesis that (a) psoriasis may reduce KRT6B methylation and upregulate its expression, thereby potentially protecting against melanoma, (b) while it may also induce melanoma development by downregulating TSPAN8 and upregulating CHI3L2 expression through increased TSPAN8 methylation and decreased CHI3L2 methylation, respectively.Two-pathway schematic linking psoriasis lesions to KRT6B, TSPAN8, CHI3L2 methylation, expression and melanoma.A) Protective Pathway: A psoriasis lesions icon on the left connects by an arrow to a hexagon labeled KRT6B hypo. A rightward arrow connects KRT6B hypo to a circle labeled KRT6B express.
A line from KRT6B express connects to a label Prevent and ends at the cutaneous melanoma icon on the right with an inhibition bar at the melanoma end. The letter a is shown near the upper left of this pathway. B) Tumor-Promoting Pathway: A psoriasis lesions icon on the left connects by separate arrows to two hexagons labeled TSPAN8 hyper and CHI3L2 hypo.
A rightward arrow connects TSPAN8 hyper to a circle labeled TSPAN8 express and another rightward arrow connects CHI3L2 hypo to a circle labeled CHI3L2 express. From each expression circle, a line connects to a label Induce and continues to the cutaneous melanoma icon on the right, ending with arrowheads at the melanoma. The letter b is shown near the lower left of this pathway.
Legend: Methylation state, Gene expression, Activation, Inhibition.
This study represents, to our knowledge, the first integrative analysis combining DNA methylation and transcriptomic data to identify cross-talk molecules between psoriasis and CM, and further provides potential diagnostic and therapeutic targets. Several limitations of this study should be acknowledged. First, because this study uses public data and the two diseases are studied based on independent datasets, causality remains elusive.
In the absence of wet-lab validation, we leveraged external transcriptomic data from TCGA and GTEx, together with IHC staining from the HPA, to support our results. Second, we did not account for CM clinical stages, which should be addressed in future work. Third, the relatively scarce data limited our external validation, which was therefore restricted to transcriptome-based CM.
Finally, potential confounding factors, such as patients’ medication use, may have influenced the results of this study.
Conclusion
In summary, this integrative transcriptomic and epigenomic analysis identified KRT6B, CHI3L2, and TSPAN8 as candidate co-regulated MRGs in psoriasis and CM. These genes suggest potential disease-specific methylation-expression patterns oppositely associated with melanoma risk—among which KRT6B and CHI3L2 merit particular attention. However, these associations are derived from public data and are insufficient to establish causality; the findings therefore require further experimental and functional validation in future studies.