{'fraction_non_alphanumeric': 0.03802698360070392, 'fraction_numerical': 0.01915376571705476, 'mean_word_length': 5.218873909595559, 'pattern_counts': {'":': 0, '<': 26, '<?xml version=': 0, '>': 2, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Cancer Associated Fibroblasts-Related Six-Gene Panel for Anti-PD-1 Therapy in Melanoma Driven by Weighted Correlation Network Analysis and Supervised Machine Learning

Background Melanoma is a highly aggressive skin cancer with a poor prognosis and mortality. Immune checkpoint blockade (ICB) therapy (e.g., anti-PD-1 therapy) has opened a new horizon in melanoma treatment, but some patients present a non-responsive state. Cancer-associated fibroblasts (CAFs) make up the majority of stromal cells in the tumor microenvironment (TME) and have an important impact on the response to immunotherapy. There is still a lack of identification of CAFs-related predictors for anti-PD-1 therapy, although the establishment of immunotherapy biomarkers is well underway. This study aims to explore the potential CAFs-related gene panel for predicting the response to anti-PD-1 therapy in melanoma patients and elucidating their potential effect on TME. Methods Three gene expression datasets from melanoma patients without anti-PD-1 treatment, in a total of 87 samples, were downloaded from Gene Expression Omnibus (GEO) as the discovery sets (GSE91061) and validation sets (GSE78220 and GSE122220). The CAFs-related module genes were identified from the discovery sets by weighted gene co-expression network analysis (WGCNA). Concurrently, we utilized differential gene analysis on the discovery set to obtain differentially expressed genes (DEGs). Then, CAFs-related key genes were screened with the intersection of CAFs-related module genes and DEGs, succeeded by supervised machine learning-based identification. As a consequence of expression analysis, gene set enrichment analysis, survival analysis, staging analysis, TME analysis, and correlation analysis, the multidimensional systematic characterizations of the key genes were uncovered. The diagnostic performance of the CAFs-related gene panel was assessed by receiver operating characteristic (ROC) curves in the validation sets. Eventually, the CAFs-related gene panel was verified by the expression from the single-cell analysis. Results The six-gene panel associated with CAFs were finally identified for predicting the response to anti-PD-1 therapy, including CDK14, SYNPO2, TCF4, GJA1, CPXM1, and TFPI. The multigene panel demonstrated excellent combined diagnostic performance with the area under the curve of ROC reaching 90.5 and 75.4% ~100% in the discovery and validation sets, respectively. Conclusion Confirmed by clinical treatment outcomes, the identified CAFs-related genes can be used as a promising biomarker panel for prediction to anti-PD-1 therapy response, which may serve as new immunotherapeutic targets to improve survival outcomes of melanoma patients.

Background: Melanoma is a highly aggressive skin cancer with a poor prognosis and mortality. Immune checkpoint blockade (ICB) therapy (e.g., anti-PD-1 therapy) has opened a new horizon in melanoma treatment, but some patients present a non-responsive state. Cancer-associated fibroblasts (CAFs) make up the majority of stromal cells in the tumor microenvironment (TME) and have an important impact on the response to immunotherapy. There is still a lack of identification of CAFs-related predictors for anti-PD-1 therapy, although the establishment of immunotherapy biomarkers is well underway. This study aims to explore the potential CAFs-related gene panel for predicting the response to anti-PD-1 therapy in melanoma patients and elucidating their potential effect on TME.
Methods: Three gene expression datasets from melanoma patients without anti-PD-1 treatment, in a total of 87 samples, were downloaded from Gene Expression Omnibus (GEO) as the discovery sets (GSE91061) and validation sets (GSE78220 and GSE122220). The CAFs-related module genes were identified from the discovery sets by weighted gene co-expression network analysis (WGCNA). Concurrently, we utilized differential gene analysis on the discovery set to obtain differentially expressed genes (DEGs). Then, CAFs-related key genes were screened with the intersection of CAFs-related module genes and DEGs, succeeded by supervised machine learning-based identification. As a consequence of expression analysis, gene set enrichment analysis, survival analysis, staging analysis, TME analysis, and correlation analysis, the multidimensional systematic characterizations of the key genes were uncovered. The diagnostic performance of the CAFs-related gene panel was assessed by receiver operating characteristic (ROC) curves in the validation sets. Eventually, the CAFs-related gene panel was verified by the expression from the single-cell analysis.

INTRODUCTION
as a routine clinical biomarker (6). Tumor mutation burden (TMB), microsatellite instability (MSI), and molecular subtypes can predict the clinical response of patients (7). Unfortunately, these detection ways are inconvenient and expensive because of complex molecular detection methods. Therefore, we necessary explore faster and more economical predictors for anti-PD-1 therapy response.
The tumor microenvironment (TME), a complex system composed of tumor cells, stromal cells, infiltrating immune cells, and abundant extracellular matrix (ECM), has a strong connection with tumor immunosuppression (8). Cancer-associated fibroblasts (CAFs) constitute vastly heterogeneous stromal cells and are prominent components of the microenvironment in solid tumors. To date, it has been recognized as one of the most promised biomarkers of immunosuppressive TME (9).
CAFs are the activated fibroblasts in the early stage of tumor development marked by the heterogeneity of origin, phenotype, and function. Studies identified CAFs promote the development of malignant tumors through various mechanisms, including secreting growth factors, remodeling the extracellular matrix, promoting angiogenesis, and mediating tumor-promoting inflammation. Furthermore, it is also involved in tumor immune escape by affecting immune cell infiltration, inducing immunosuppression, and inhibiting lymphocyte tumorkilling effects (10,11). Studies demonstrated that CAFs can inhibit CD8 + T cells and activate FoxP3 + lymphocytes by secreting IL-6 and thus promote immunosuppression of TME (9). The tumor-associated macrophages (TAM) are part of inflammatory infiltrating cells in TME, of which the M2 type has immunosuppressive functions. CAFs have been proved to recruit monocytes and induce their differentiation into M2 type via monocyte chemotactic protein 1 (MCP-1) and chemokine CXCL-12 (12). Importantly, CAFs have been proved to directly reduce the activation of CD8 + T cells and natural killer (NK) cells by upregulating the expression of immune checkpoint signals (e.g., PD-1, PD-L1, PD-L2), which in turn hinders anti-tumor immunity (13,14). In conclusion, CAFs can promote tumor immunosuppression and immune escape by interacting with the immune cell through the secretion of cytokines, chemokines, and other cytokines.
The evidence above suggests the association between CAFrelated genes and the effectiveness of anti-PD-1 therapy in tumor patients. There are some studies foucusing on biomarkers for the response prediction of anti-PD-1 therapy in melanoma, but the majority of them tend to ignore the exploration of the functional correlation between genes as well as the regulatory mechanism of genes in TME. Furthermore, the exploration of anti-PD-1 therapy response biomarkers related to CAFs in melanoma is few up to now. Therefore, in this study, a robust analysis, comprised of weighted gene co-expression network, least absolute shrinkage and selection operator (LASSO), and random forest machine learning, was applied to screen CAFs-related gene panel associated with anti-PD-1 therapeutic response in melanoma, aiming to understand the potential functions of these genes in TME and evaluate their anti-PD-1 therapeutic response prediction performance.
Here, three existing RNA-seq datasets on anti-PD-1 immunotherapy of melanoma were analyzed using a bioinformatics approach. First, GSE91061 was applied to construct a gene co-expression network by weighted gene co-expression network analysis (WGCNA) to screen key modules related to CAFs score. Functional enrichment analysis was further to investigate the biological function of the key module genes. Second, six CAFs-related genes for predicting response to anti-PD-1 therapy were identified with the combination of LASSO and random forest. Subsequently, adopting gene set enrichment analysis (GSEA) and gene expression profiling interactive analysis (GEPIA), the biological functions and survival/staging associations of the six genes were explored based on TCGA-SKCM RNA sequencing data. To gain insight into the impact and association of these molecules on immune/stromal cells in TME, xCELL was applied to explore the correlation of the six genes with immune/stromal cells. Not only that, the multivariate Cox model confirmed the excellent performance of the six-gene panel in survival prediction, and TME analysis demonstrated significant TME alterations existed in different risk groups. Further gene correlation analysis revealed significant positive correlations of these six genes with anti-tumor immune-related molecules, immune checkpoint-associated molecules, and CAFs-related stimulators, which suggested the potential anti-tumor immune functions. After validating the prediction accuracy of the gene panel in the GSE78220 and GSE122220 datasets, we finally summarized key features of the KCDEGs and investigated single-cell localizations to determine their significant expression in CAFs and TME-associated cells.

METHODS
The workflow of this study is shown in Supplementary Figure 1, including three main parts, WGCNA-based CAFs-related differential gene screening, key gene extraction by machine learning models, and key gene validation based on multidimensional system analysis. The details are described as follows.

Construction of Weighted Co-expression Network in Melanoma Patients
WGCNA aims to mine functionally related genes with similar co-expression patterns (22). Expression correlation coefficients were applied to measure intergenic co-expression relationships, and genes with high correlation coefficients were assigned to the same module and presented similar expression patterns. These highly correlated module genes may be involved in the same biological process or pathway. Gene co-expression networks were constructed using the WGCNA package in R (4.1.1) by the top 25% of the median absolute deviation (MAD) ranked genes in the GSE91061 expression matrix normalized by log2(count+1). The co-expression similarity matrix for all genes was generated using the average linkage method and subsequently converted into an adjacency matrix to ensure the construction of an unsigned scale-free network. The weighted adjacency matrix was further transformed into a topological overlap matrix (TOM) to estimate the connectivity of the unsigned network. A hierarchical clustering approach was adopted to construct the clustering tree structure of the TOM, and dynamic hybrid cuts were achieved by the variability of the TOM with a cut height of 0.25. Gene modules were fused under the variability of the estimated module eigengenes (MEs) and represented by different colors.  (14). Then, according to the marker genes, the CAFs score was calculated for each patient in GSE91061 (n = 49) using the GSVA algorithm. Finally, the correlation between the CAFs score and MEs was assessed using Pearson's correlation coefficient, and the gene module with the smallest P-value was selected as the most relevant module for the CAFs score. Genes in the module were defined as highly correlated with CAFs.

Identification of Key Genes by LASSO Regression Analysis and Random Forest Algorithm
RNA-seq count data from GSE91061 were normalized with TMM (trimmed mean of M values) and analyzed using the edgeR package (28), and genes with P < 0.05 were defined as anti-PD-1 response differentially expressed genes (DEGs). The intersection between module genes associated with CAFs and DEGs related to anti-PD-1 therapy response were extracted and defined as CAFsrelated DEGs (CDEGs), and further LASSO and random forest analysis were adopted to identify key CDEGs (KCDEGs) using the glmnet package (nlambda = 1000, 10-fold cross-validation) and randomForest package. LASSO regression was characterized by variable selection and regularization while fitting a generalized linear model, which aids it well-adapted to the linear and nonlinear operations (29) and random forest is a popular classifier that can handle input samples without dimensionality reduction to generate an unbiased estimate of the error in the process of forest building (30).

Gene Set Enrichment Analysis for KCDEGs
Once 458 TCGA-SKCM samples were divided into two groups with the median gene expression, the edgeR package was further employed for differentially expressed genes (FDR < 0.05) between the two groups. Subsequently, based on hallmark gene sets of the MSigDB (http://www.gsea-msigdb.org/gsea/msigdb/ index.jsp) (31), GSEA was utilized (FDR < 0.05) to estimate the function of KCDEGs by clusterProfiler package (32).

Survival and Staging Analysis of KCDEGs in Melanoma Patients
To investigate whether KCDEGs have an impact on survival and disease progression in melanoma patients, we applied gene expression analysis (One-way ANOVA, P < 0.05), survival analysis (Log-rank P < 0.05), and staging analysis [One-way ANOVA, Pr(>F) < 0.05] in TCGA-SKCM patients by GEPIA (33). In addition, we evaluated the predictive effect of combined KCDEGs on patients' prognosis in the TCGA-SKCM cohort (n = 458). A multivariate Cox proportional risk regression model was constructed for the KCDEGs using the survival R package (34), as well as regression coefficients were retained for each KCDEG. And risk score was calculated for each patient according to the following equation.
SKCM patients were assigned to low-and high-risk groups based on the median risk score calculated above. We finally assessed the overall survival between the two risk groups by the Kaplan-Meier curve.

Association Assessment of KCDEGs With TME Signature
The RNA-seq count matrix of TCGA-SKCM was converted into TPM (transcripts per million) values for the next analysis. The variation of TME for each sample and 64 xCell TME signature were investigated by the xCell online tool (https:// xcell.ucsf.edu/) (35) to estimate the infiltration level of TMEassociated cells (mainly include immune cells and stromal cells). Pearson correlation analysis was utilized to calculate the degree of correlation between KCDEGs and the infiltration levels of TME-associated cells (P < 0.05). We used the Wilcoxon test to compare the TME-associated cell infiltration level between lowand high-risk groups (P < 0.05).

Correlation Analysis of CAFs-Related Stimulators and Immunotherapy-Related Genes
The mRNA expression data for KCDEGs, CAFs-related stimulators, antitumor immune-related genes and, immune checkpoint-related genes were extracted from the TPM matrix of TCGA-SKCM. Spearman correlation analysis was adopted to calculate the correlation coefficient between each gene (P < 0.05). And the ggcorrplot R package (36) was employed for the correlation coefficient visualization.

Diagnostic Performance Evaluation and Validation of KCDEGs
To further assess the predictive ability of KCDEGs for anti-PD-1 treatment response, the pROC package (37) was chosen to calculate the area under the curve (AUC) of the receiver operating characteristic (ROC) curve for each KCDEGs in discovery and validation sets. Higher AUC values represent a better gene diagnostic performance. By binary logistic regression (SPSS, v25.0), KCDEGs were firstly fitted to the CAF-derived gene panel, and the diagnostic performance of the CAF-derived gene panel was assessed through ROC curves. To determine the single-cell expression distribution of these KCDEGs, we analyzed the expression of each KCDEGs with single-cell RNA-seq data by Tumor Immune Single Cell Hub (TISCH) (38). Six SKCM datasets were included in the analysis: four non-treatment cohorts (GSE123139, GSE139249, GSE148190 and GSE72056), one anti-PD-1 treatment cohort (GSE115978) and one anti-PD-1/anti-CTLA-4 treatment cohort (GSE120575).

Statistical Methods
Statistical analyses and visualization were performed with R software (version 4.1.1). Normally distributed continuous variables between two groups were compared by t-test. Otherwise, the Wilcoxon test was applied. The P-values for GSEA were corrected by the Benjamini-Hochberg method. The statistical significance of survival analysis was assessed by the Log-rank test and corrected by the Bonferroni method. Correlations between variables were explored by Pearson or Spearman coefficients. All statistical tests were two-sided.

Construction of Weighted Gene Co-expression Networks and Identification of CAFs-Related Modules
To obtain the gene modules associated with CAFs, we constructed a weighted gene co-expression network by WGCNA for GSE91061. A soft threshold of β = 10 ( Figure 1A) was chosen for the scale-free topological network model fitting (R2 = 0.85). Subsequently, 14 gene modules were obtained by average linkage hierarchical clustering, TOM dynamic hybrid cut (cutHeight = demonstrated the similarities and differences among these modules. Similar module genes may be involved in similar biological function regulation. Finally, the Pearson correlation analysis was performed to calculate the correlation of MEs with the CAFs score and clinical characteristics (Non-response and Response). The results suggest that the green module was the most relevant to the CAFs score, which included 184 module genes ( Figure 1E). Further, we displayed a significant correlation between module membership in the green module and gene significance for CAFs ( Figure 1F, Cor = 0.89, P = 5.6E-64).

Biological Functions of Genes in the Green Module
The green module was identified as the most associated gene module with CAFs. To further explore the functional association between green module genes and CAFs, enrichment analysis (FDR < 0.05) of module genes was realized based on GO, KEGG, REAC, and WP. GO results exhibited that the molecular function of 184 genes is mainly enriched in extracellular matrix structural constituent, structural molecule activity, and growth factor binding (Figure 2A). The involved terms of biological processes and cellular components include extracellular matrix organization, extracellular structure organization, extracellular matrix, and collagen-containing extracellular matrix. The above results demonstrated the significant association of module genes with ECM remodeling. Increasing evidence has revealed that CAFs are major contributors to ECM remodeling in the TME (14,39,40). CAFs can alter the structure, arrangement, and stiffness of the ECM by secreting various matrix proteins, growth factors, and cytokines, which in turn affect the migration and invasion of cancer cells (41). Another study emphasized that CAFs can also promote the expression of fibronectin and laminin through the secretion of cytokine TGF-β1, which affects the remodeling process of ECM (42). Subsequent pathway enrichment results of KEGG, REAC, and WP similarly demonstrated the close association of these module genes with ECM (Figures 2B-D), such as ECM-receptor interaction, ECM proteoglycans, and miRNA targets in ECM, and membrane receptors. Notably, the enrichment results for KEGG and WP presented significant associations of these genes with the PI3K-AKT signaling pathway (KEGG: FDR = 2.4E-05; WP: FDR = 4.9E-04), suggesting that CAFs-related genes may be involved in the regulation of the PI3K-AKT signaling pathway. The previous study has revealed that CAFs can regulate proliferation, apoptosis and, invasion of lung cancer cells by activating the PI3K-AKT-mTOR signaling pathway through secretion of IL-22 (43).

Identification of Key CAFs-Related
Differential Genes by Machine Learning 2191 DEGs were retained by performing gene differential analysis of anti-PD-1 treatment-responsive and non-responsive patients in the GSE91061 dataset. Subsequently, we obtained 27 CDEGs with the intersection of 184 CAFs-related genes and 2191 DEGs (Figure 3A). To further identify the KCDEGs, a linear model was first constructed for the 27 CDEGs using LASSO regression with variable screening and regularization ( Figure 3B). The meansquared error of the model was minimized when log(λ) = −2.788 ( Figure 3C). Finally, we screened nine genes as the best variables to distinguish patients responding and non-responding to anti-PD-1 treatment. Concurrently, a random forest classifier was built based on 27 CDEGs, and we selected ntree = 150 to ensure the stability of the constructed model ( Figure 3D). Then all molecules were ranked by mean decrease accuracy and mean decrease Gini (Figure 3E). The intersection of the top20 genes ranked by two methods in the random forest model and the 9 genes identified by LASSO was extracted. Finally, we obtained 6 KCDEGs identified by two machine learning algorithms ( Figure 3F).

Expression Analysis of KCDEGs
KCDEGs may display diversified expression patterns in different patients, and we further investigated their expression between patients responding and non-responding to anti-PD-1 therapy. Supplementary Figure 2A showed that CDK14 (Cyclin Dependent Kinase 14), SYNPO2 (Synaptopodin 2), TCF4 (Transcription Factor 4), and TFPI (Tissue Factor Pathway Inhibitor) expression were significantly elevated in the response group (Student's t-test, P < 0.05). Although the expression of GJA1 (Gap Junction Protein Alpha 1) and CPXM1 (Carboxypeptidase X, M14 Family Member 1) had no significant change, CPXM1 tended to increase somewhat in the non-response group (Student's t-test, P = 0.054). In addition, we investigated the expression of KCDEGs in GSE114445 (Supplementary Figure 2B). Among them, expressions of CDK14, SYNPO2, TCF4, and GJA1 displayed a gradual decrease in normal skin, nevus and melanoma samples. In contrast, the expression of TFPI was significantly increased in melanoma samples compared to nevus (Student's t-test, P < 0.05), but without change in normal skin (Student's t-test, P = 0.73). Notably, CPXM1 expression demonstrated a progressive or gradual increasing trend in normal skin, nevus, and melanoma samples while decreasing in patients responding to anti-PD-1 therapy, which showed the gene might serve as a pro-cancer molecule in melanoma.

Gene Set Enrichment Analysis of KCDEGs
As shown in Figure 4A, all enriched KCDEGs in 28 cancer hallmarks (FDR < 0.05) were involved in their respective biological processes in a certain state. For example, the enriched KCDEGs-related genes involved in DNA repair, MYC targets, E2F targets, and oxidative phosphorylation were repressed, while others were activated in epithelialmesenchymal transition (EMT), KRAS signaling up, allograft rejection, and inflammatory response. Moreover, EMT was the most significant activated term and closely associated with all KCDEGs ( (Figure 4B). Taken together, the above results fully demonstrate the close association of KCDEGs with cancer development, especially the regulatory association with EMT confirmed in several studies (44)(45)(46). In addition, previous studies have identified an association between EMT and elevated PD-1 levels in lung adenocarcinoma patients (47). This evidence suggests a potential link among CAFs, EMT, and anti-PD-1 therapy response. The Role of KCDEGs in the Melanoma TME TME was suggested to have an important effect on the progression of malignant tumors and response to anti-PD-1 therapy (48,49), which is confirmed by our findings. Twenty-one TME signature fractions significantly correlated with KCDEGs were traced in our further investigation (Pearson correlation, P < 0.05). As demonstrated in Figure 6A, KCDEGs were significantly and positively correlated with stromal fractions such as fibroblasts, endothelial cells, and mesangial cells, but a significant negative correlation with melanocytes. In the analysis of the relationship between KCDEGs and immune fractions, TCF4 and TFPI were positively associated with immune score  (P < 0.01), while CPXM1 showed a negative association (P < 0.05). Meanwhile, TCF4 and TFPI were positively related to CD4 + memory T cells, dendritic cells (DCs), monocytes, CD4 + T cells (P < 0.05) and negatively related to Th1 and macrophages M2 (P < 0.05). CPXM1 was only positively correlated with CD4 + memory T cells (P < 0.05) and negatively correlated with macrophages (P < 0.05).

Survival Impact and Staging Expression of KCDEGs on Melanoma Patients
The complexity of the relationship between KCDEGs and TME signature also indicates that they act as different mechanisms in TME. To better evaluate the ability of combined KCDEGs impacting on survival and TME, we performed multivariate Cox regression analysis and calculated risk score according to the following equation: +0.095 * Exp CPXM1 + 0.072 * Exp TFPI Kaplan-Meier curve demonstrates that the high risk score group (high-risk group) has worse survival (Log-rank P = 0.00011, Cox P = 6E−05, HR = 2.7, 95% CI 1.7-4.4) (Figure 6B). In addition, there is a significant difference in TME between the high-and low-risk groups (Figures 6C,D). In particular, compared with the low-risk group, the components of the immune score and macrophages, monocytes, as well as CD8 + naïve T cells are lower (Wilcoxon test, P < 0.0001), while the components of melanocytes, Tregs, and CD4 + memory T cells are higher in the high-risk group (Wilcoxon test, P < 0.0001).
We can conclude the following: (i) KCDEGs may affect the secretion of CAFs-related stimulating factors or participate in the activation of CAFs. For example, IL-1β and IL-6 are proinflammatory cytokines that can promote the activation of CAFs through NF-κB and JAK-ROCK-STAT3 signaling pathways, thus causing tumor progression (51,52). Heat shock factor 1 (HSF1) regulates the pro-tumor effects of CAFs by activating β-catenin and YAP/TAZ signaling pathways (53). (ii) Antitumor immunerelated molecules including CXCL10, CXCL9, CXCR3, CCL5,  Heatmap and box plot depicting differences of TME features in two risk score groups. Wilcoxon test, *, P < 0.05; **, P < 0.01; ***, P < 0.001; ****, P < 0.0001. that they may affect the benefit from ICB therapy and turn into potential targets for immunotherapy.

Evaluation and Validation of Response Prediction Performance of KCDEGs for Anti-PD-1 Treatment
The above results of the correlation analysis have suggested a strong association between KCDEGs and ICB therapy. Six KCDEGs (except SYNPO2, whose corresponding probe was not found in GSE122220) were evaluated in discovery and validation sets to predict the response to anti-PD-1 therapy by ROC curves. The predictive power of the ROC was inferred by calculating its AUC, with a larger AUC indicating better prediction. In the discovery set, CDK14, SYNPO2, TCF4 and TFPI are the best predictors (AUC > 70%) (Figure 7A). In the two validation sets, the AUC ranges from 56.4 (SYNPO2 in GSE78220) to 75% (TCF4 in GSE12220). Overall,
The results of the single-cell analysis demonstrated that CDK14 is predominantly distributed in B cells (Figure 8B), SYNPO2 presented highest expression in fibroblasts (Figure 8C), TCF4 mainly enriched in fibroblasts, B cells, and endothelial cells (Figure 8D), GJA1 presented highest expression in fibroblasts and endothelial cells (Figure 8E), CPXM1 mainly distributed in fibroblasts (Figure 8F) while TFPI exhibited the highest expression in fibroblasts and endothelial cells (Figure 8G). Additionally, the results were verified in the other five datasets (Supplementary Figure 4).

DISCUSSION
Although anti-PD-1 therapy is expensive and exerts adverse effects, anti-PD-1 immune checkpoint inhibitors have been confirmed of excellent clinical activity in melanoma treatment (55). However, not all patients can produce an objective response to it (4). It is essential and urgent to determine the patients' response using reliable biomarkers to predict anti-PD-1 efficacy before the treatment. CAFs can suppress antitumor immune responses through upregulation of PD-1 expression, and CAFsrelated genes may be the potential markers in predicting the efficacy of anti-PD-1 immune checkpoint inhibitors (56)(57)(58). In this study, we successfully elucidated the potential regulatory efficacy of CAFs-related genes in TME and discovered a sixgene panel for prediction of response to anti-PD-1 therapy. Our discovery manifests the panel is favorably consistent with the clinical response of patients, which can aid in the better clinical management of melanoma patients.
We first calculated the CAFs score for each patient based on 56 CAFs markers reported in literature through the GSVA algorithm. Subsequently, WGCNA was applied to construct gene co-expression networks, which has the potential to identify CAFs-related genes (59),and the green module from networks was of the highest correlation with the CAFs score. By intersecting green module genes and the DEGs, we obtained 27 response-associated CDEGs for anti-PD-1 therapy. Finally, 6 KCDEGs (CDK14, SYNPO2, TCF4, GJA1, CPXM1, and TFPI) were identified by the methods of lasso regression and random forest.
Comprehensive analysis implied that KCDEGs are closely associated with CAFs and appeared as potential predictors of anti-PD-1 therapy efficacy in melanoma. First, gene expression analysis demonstrated six KCDEGs expressed with significant differences could reflect patients' response status to anti-PD-1 therapy, with CDK14, SYNPO2, TCF4, TFPI up-regulated in responders. In addition, the expression of six KCDEGs presented significant differences in SKCM. Secondly, GSEA indicated that KCDEGs were enriched in several cancer hallmarks such as EMT, MYC targets, KRAS signaling, and inflammatory response. EMT is the process epithelial cells acquire a mesenchymal phenotype after downregulating epithelial features, which exhibited crucial functions in cancer progression. Previous studies confirmed that EMT could promote CAF activation and alter the secretory phenotype of CAFs in the mesenchyme (60)(61)(62). Additionally, MYC is one of the widely studied oncogenes, and high MYC targets scores in breast cancer were associated with increased tumor immune cell infiltration (63). Moreover, CAFs have been verified to have complex crosstalk in KRAS-driven tumorigenesis and to influence immune regulation through the immune checkpoint (64). Taken together, KCDEGs may be a vital part of the melanoma immune infiltration regulatory network and suggest a predictive effect on the response to anti-PD-1 therapy.
Survival analysis shows that high expression of CDK14 and SYNPO2, plus low expression of GJA1 and CPXM1 were associated with better survival in SKCM patients. As patients with excellent prognoses often displayed positive responses to anti-PD-1 therapy (65), we speculated that the combination of CDK14 and SYNPO2 (high expression), as well as GJA1 and CPXM1 (low expression), may be associated with a positive response to anti-PD-1 therapy. As an activator of the WNT pathway through mediating the phosphorylation of LRP6 in the G2/M phase, CDK14 was found that it could promote the progression of several cancers (66)(67)(68). However, our findings suggest that high expression of CDK14 is associated with better prognosis and anti-PD-1 therapy response in melanoma patients, implying a unique mechanism of CDK14 in melanoma. As expected, SYNPO2, whose promotor methylation and low transcriptional expression were found to associate with metastasis and poor clinical outcome in melanoma (69), can be a potential predictor to provide new therapeutic strategies in this study. GJA1 is a gene encoding the protein connexin 43 (Cx43) that is a vital component of gap junctions and was disclosed to be of an important role in communication between tumor cells and surrounding immune cells like NK and dendritic cells (DC) (70). The up-regulation of GJA1 expression implicates a crucial link between melanoma cells and endothelial cells, thus enhancing tumor metastasis in melanoma (71). Interestingly, Cx43 has generally been reported as a potential cancer suppressor to improve therapeutic efficacy against melanoma (72). CPXM1 is suggested as an immune-related gene to predict cancer prognosis (73,74), which is consistent with our findings.
Correlation analysis of KCDEGs with various fractions in TME revealed all these six KCDEGs (TFPI, TCF4, CDK14, SYNPO2, GJA1, and CPXM1) positively correlated with stroma score, endothelial cells, and fibroblasts ( Figure 6A), which can be promisingly accountable for tumor progress and prognosis of melanoma patients. For example, TFPI, encoding tissue factor pathway inhibitor, is closely associated with the endothelial cell surface and can regulate angiogenesis (75). Notably, endothelial cell activation requires the elevated expression of Snail1, which can strengthen CAF activation in a paracrine manner promoting immune infiltration in TME by secreting IL-6 in tumor development (9,76). All of these are witnessed by the coculture of CAFs with endothelial cells under hypoxic conditions promoting breast cancer angiogenesis (77) and tumor exposedlymphatic endothelial cells (teLEC) reported as the promoter of cancer cell invasion and tumor cell proliferation by regulating IL-6 (78). Taken together, we believe TFPI may mediate the action of CAFs and endothelial cells through transcription factors and cytokines, thereby regulating immune cell infiltration in TME, which is also investigated from the correlation analysis of KCDEGs-immune cell fractions. TCF4 and TFPI, up-regulated in the PD-1 therapy response group and low expressed in the SKCM group (Figures 5A,B), were positively related to CD4+ memory T cells, CD4+ T cells, monocytes, DC cells ( Figure 6A). In other words, the decreased expression of TCF4 and TFPI are associated with a low level of immune infiltrations or immune desert in TME. Interestingly, current studies keep the point that a certain of the desert immune microenvironment (DIM) in patients caused by CAFs can result in poor outcomes of anti-PD-1 therapy for patients (11). It is reasonable that TCF4 and TFPI play an important role in inhibiting CAFs-mediated DIM. For instance, TCF4, known as immunoglobulin transcription factor 2 (ITF-2; E2-2; SEF-2), contributes to the development of lymphoid (79) and mature plasmacytoid dendritic cells (pDCs), including Pitt-Hopkins syndrome caused by E2-2 haploinsufficiency (80) and the tumor-suppressive effect in SHH medulloblastoma (81). In all, the inhibition of TCF4 and TFPI expression may lead to CAFs-related desert immune microenvironment, which results in the non-response to anti-PD-1 therapy of patients.
We also analyzed the correlation between six KCDEGs and CAFs-related stimulators, antitumor immune-related molecules as well as immune checkpoint-related molecules. All KCDEGs except CPXM1 and GJA1 were positively relevant to most of the mentioned molecules. This suggests that KCDEGs are directly or indirectly involved in the activation process of CAFs and affect the response to immunotherapy in melanoma patients, which is consistent with our previous conclusion that CAFsrelated genes can serve as a predictive marker panel of response to anti-PD-1 therapy.
Although the six KCDEGs can distinguish responders and non-responders of anti-PD-1 immunotherapy, the therapeutic efficacy is often influenced by multiple factors in the TME and immune system. Our results also demonstrated an inconsistent correlation of each KCDEGs with immune fractions. This implies the predictive limitation of the accuracy and specificity in anti-PD-1 immunotherapy using a single marker. Combining two or more markers may be more effective for predicting the efficacy of anti-PD-1 therapy (82). To this end, we first employed multivariate Cox regression analysis to assess the OS and TME variety of combined-KCDEGs in SKCM patients. Survival analysis revealed worse OS in the high-risk KCDEGs group. In addition, the two groups had significant differences in multiple TME metrics and anti-tumor immune components. Immune score and cells (including CD8+ naive T cells, DCs, macrophages, and monocytes) tended to decrease in the highrisk group, suggesting that the overall high-risk score of KCDEGs relates to immune infiltration suppression and PD-1 therapy non-response. And this is supported by poorer survival in the high-risk group. Surprisingly, the AUC of combined-KCDEGs (the six-gene panel) was 90.5% in the discovery sets and reached 75.4∼100% in the validation sets, which is a substantial improvement compared with the previous study (83) (AUC = 0.75 in the discovery sets and 0.71 in the validation sets). These results fully demonstrate the significant advantage of KCDEGs in anti-PD-1 efficacy prediction. It also indicates that the vital role of CAFs in patient response to anti-tumor immunotherapy may be beyond our imagination and deserves to be explored in detail.
We identified six key genes with predictive value for anti-PD-1 treatment efficacy in melanoma patients via WGCNA and supervised machine learning. However, several limitations remain in this study. We applied robust computational biology methods, including WGCNA, supervised random forest, and LASSO regression, to screen predictive gene panels from multiple perspectives, but validation by experimental biology methods such as RT-PCR is also urgently needed. Secondly, the limited number of samples may make the efficacy assessment model inaccurate, cohorts with small samples often lead to overfitting of machine learning model and reduce the prediction accuracy, thus the model needs to be validated in a larger patient cohort. Thirdly, it is necessary to evaluate the actual predictive effect of the sixgene panel in a clinical patient cohort. Finally, the regulatory mechanism of KCDEGs induced by CAFs in immunotherapeutic suppression is unclear and deserves to explore.

CONCLUSION
Using WGCNA combined with supervised machine learning algorithms, we identified a novel CAFs-related panel, including six genes (CDK14, SYNPO2, TCF4, GJA1, CPXM1, TFPI), which can distinguish the response of melanoma patients under anti-PD-1 immunotherapy. The multigene may become a potential biomarker panel to guide immunotherapy in the future.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
LT and FL contributed to the study inception and design. LT, FL, JL, and YT equally analyzed the data. JC and ML contributed to the study design and study supervision. All authors contributed to the manuscript writing and approved the final version of the manuscript.

ACKNOWLEDGMENTS
We would like to thank the peer reviewers input in improving this manuscript.