A set of melanoma immune escape related genes were extracted from the Lawson et al. (26) ( Supplementary Material 1 ). Immune-related genes were extracted from the IMMPORT database ( Supplementary Material 1 ). Gene expression profiles, clinical follow-up information, and somatic mutation data were extracted from The Cancer Genome Consortium (TCGA) database through the ‘TCGAbiolinks’ package. In addition, we extracted normal control samples from the Genotype-Tissue Expression Project to compensate for the absence of normal sample controls in the TCGA cohort. Single cell transcriptome data as well as external validation datasets were downloaded from the Gene Expression Omnibus (GEO) and International Cancer Genome Consortium (ICGC) databases ( Table 1 ). Duplicate genes are processed by the avereps function of the ‘limma’ package. Processing of gene expression values into Transcripts Per Kilobase of exon model per Million mapped reads (TPM) and normalized by log2. With different batches of GEO datasets we remove the batch effect by using the combat function of the ‘sva’ package.

Somatic mutation data in mutation annotation format (MAF) were visualized through ‘maftools’ R package, which provides a large number of analysis and visualization modules commonly used in cancer genomic studies (27).

cBioPortal provides a Web resource for exploring, visualizing, and analyzing multidimensional cancer genomics data (28). Differentially altered proteins (both unphosphorylated and phosphorylated) associated with immune escape related genes were identified in cBioPortal by the Reverse Phase Protein Array (RPPA) module.

Protein-protein interaction network of protein-coding genes constructed by STRING database. Minimum required interaction score 0.4 and disconnected nodes in the network were hide.

The melanoma patients from TCGA cohort were divided into a training set (TRS) and a testing set (TES). The TRS was used to construct a prognostic risk model of melanoma and the TES was utilized to valid the predictive capability of this model. ICGC-SKCM-US is used as external dataset validation. Risk prognostic models were constructed by univariate cox regression, least absolute shrinkage and selection operator (LASSO) regression analysis, and multivariate cox regression analysis. The risk score for each sample is calculated as follows:

Exploring the diagnostic value of risk scores for 1-year, 3-year, 5-year, 8-year and 10-year survival status with the ‘pROC’ package.

‘GSVA’ package was used to explore the signaling pathways in which the screened protein-coding genes are involved and the c2.cp.kegg.v7.5.1.entrez.gmt gene sets was used to calculate differences in enrichment scores for pathways in different populations in the tow group.

To explore the role of our screened protein-coding genes and the tumor microenvironment, we used ESTIMATE, EPIC, TIMER, quanTiseq and MCPcounter algorithms to calculate the proportion of various immune factors infiltrating the tumor microenvironment and explored the correlation between protein-coding gene expression and immune factor infiltration levels.

Considering that the number of cells in one single cell dataset is too small, we removed the batch effect and merged GSE72056 and GSE186344 via the ‘harmony’ package. The integrated single cell data was then analyzed using the ‘Seurat’ package, including finding highly variable genes, centralization, PCA downscaling, cell clustering, tSNE (t-Distributed Stochastic Neighbor Embedding) and UMAP (Uniform Manifold Approximation and Projection) nonlinear downscaling, finding differential genes, and cell annotation. To investigate cell-to-cell interactions and to determine the mechanisms of communication molecules at single-cell resolution, 8 cell subgroups were studied using the R package ‘CellChat’.

To explore the potential of the screened protein-coding genes as predictive biomarkers for chemotherapy or immunotherapy, we attempted to assess the correlation between different expression populations of protein-coding genes and responsiveness. ‘oncoPredict’ is an R package for predicting in vivo or cancer patient drug response and biomarkers from cell line screening data (29). We extracted immunotherapy cohort (PRJEB23709 and PRJEB25780) from Tumor Immune Dysfunction and Exclusion (TIDE) to explore the potential of protein-coding genes as predictors of immunotherapy response ( Supplementary Material 2 ).

The PIG1 cells (Otwo Biotech, ShenZhen Inc. China), A2058 and SKMEL28 cells (the Chinese Academy of Sciences) were cultured in the Dulbecco’s modified Eagle’s medium (DMEM) (Gibco, Thermo Fisher Scientific, Inc.) medium containing 10% FBS (Gibco; Thermo Fisher Scientific, Inc.) in an incubator at 37˚C, and the air of the incubator consisted of 5% CO2.

The cells were cultivated in a 6-well plate at a density of 40×10⁴ cells per well and incubated at 37˚C. Following collected the cells after the cell’s density reached 80%. Total RNA was extracted from the 6-well plate using TRIzol reagent (Thermo Fisher Scientific, Inc.), then subjected to reverse transcription-quantitative polymerase chain reaction (RT-qPCR) to detect the mRNA expression of LCK. qPCR was performed with a SYBR Green Real Time PCR kit (Thermo Fisher Scientific, Inc.) on CFX96 Touch Real Time PCR System (BioRad Laboratories, Inc.). The primers used for real-time PCR were at Table 3 ( Table 2 ). qPCR was performed under the following conditions: 3 min at 95˚C for 1 cycle, 10 sec at 95˚C, 30 sec at 65˚C for 39 cycles, and 95˚C for 5 sec. Changes in the expression of target genes were calculated using the 2^-ΔΔCq method.

A set of immune-escape related genes and a set of immune-related genes were extracted from the previous study and from the IMMPORT website respectively, which contains 182 genes and 2483 genes respectively in total. Overlapping these two gene sets, 31 core immune-escape related genes ( Table 3 ) were determined as the subjects for further analysis ( Figure 1A ). In addition, 472 (including 471 tumor samples and 1 normal sample) and 398 normal samples were obtained from the TCGA and GTEx databases, respectively. Prognostic Analysis showed that except CALR, TNFRSF1A, HDAC1, JAK1 and TFRC, which were not significantly different, and MAPK1, which was an unfavorable prognostic factor, the rest of all the core immune-escape genes were favorable prognostic factors, and all of them were significantly different ( Figure 1B ). Expression analysis showed that the expression of all core immune escape related genes was dysregulated. IFNGR1, JAK2, SOCS1, IKBKG, JAK1, TNFAIP3, TNFRSF1A, FAS, IKBKG, TBK1 and TGFBR2 were highly expressed in normal tissues, while the expression of the remaining genes were upregulated in tumor tissues, all of the above differential expression analysis results were significantly different (p-value< 0.05) ( Figure 1C ). Genes with significantly different in both expression difference analysis and prognostic analysis were initially screened out for inclusion in the next analysis. Somatic mutation profiles of 467 melanoma patients downloaded from TCGA database were analyzed and visualized via the ‘maftools’ R package (30). The results showed that they all had low mutation rates ( Figure 1D ). In addition, identification of differentially expressed proteins regulated by primary screening genes through the RPPA module of the cBioPortal database and 164 differentially expressed proteins were obtained in total ( Figure 2 ).

Through STRING database, protein-protein interaction network was constructed based on the above differentially expressed proteins, including 188 nodes and 3670 edges in summary. Subsequently, we filtered out the top 60 nodes in the entire network in terms of connectivity ( Figure 3A ). Meanwhile, the melanoma samples were divided into TRS and TES in a roughly 1:1 ratio (228 samples in TRS and 226 samples in TES). The result of univariate cox regression analysis of TRS showed that a total of 15 prognostic factors were determined. Except CDKN1B, LCK, and RICTOR which were favorable prognostic factors, the rest of these prognostic factors were associated with reduced overall survival (p-value< 0.05) ( Figure 3B ). Immediately, LASSO regression analysis of prognostic factors was performed and 10 representative protein-coding genes were identified ( Figures 3C, D , Table 4 ). Multivariate cox regression analysis was conducted within these 10 representative protein-coding genes and 4 independent prognostic factors were finally identified which were related to prognosis in melanoma ( Figure 3E ). Overall, all 4 independent prognostic factors were associated with reduced overall survival, except for LCK, which was favorable factor ( Figure 3E ). Overlapping the top 60 nodes in the protein-protein interaction network in terms of connectivity and the 4 significant independent prognostic factors, we obtained 3 key prognostic factors, respectively KIT, EGFR and LCK ( Figure 3F ).

Through analyzing the independent prognostic factors in TRS, a risk prognostic model was constructed in TRS, using which the prognosis of patients could be effectively predicted. To verify the predictive efficacy of the model, the same model was constructed in TES and the ICGC cohort. The results of the risk factor analysis showed that the number of patient deaths clustered significantly as the risk score increased ( Figure 4A ). Survival analysis in the high- and low-risk score groups showed that the high-risk group was associated with lower overall survival (p-value< 0.05) ( Figure 4B ). Finally, time-independent ROC curves were established to verify the accuracy of the model. The AUC values of the ROC curves for the three cohorts indicate that the model has good accuracy ( Figure 4C ).

The infiltration level of immune cell was calculated for each sample of melanoma using four methods, EPIC, TIMER, quanTiseq and MCPcounter, respectively. Focusing on the three key independent prognostic factors mentioned above, we found that samples with high LCK expression were accompanied by higher levels of immune infiltration ( Figure 5A ). Therefore, we suggest that LCK may affect prognosis by altering the tumor microenvironment of melanoma and thereby. Survival analysis in both the TCGA cohort and the GSE54467 cohort of melanoma showed that the group with high LCK expression levels was associated with prolonged overall survival (p-value< 0.05) ( Figure 5B ). In normal and tumor tissues, differential expression analysis based on TCGA cohort and GTEX cohort showed that LCK was highly expressed in tumor tissues (p-value< 0.05) ( Figure 5C ). To verify the expression pattern of LCK, RT-qPCR was used to detect the mRNA expression of LCK. The results showed that the expression level of LCK was significantly higher in melanoma cell lines (A2058, SKMEL28) than in normal skin cell lines (PIG1, p-value< 0.05) ( Figure 5D ). We then analyzed the correlation between LCK expression levels and four clinical parameters, including tumor status, metastasis, pathological stage, and the extent and size of the primary tumor. The results indicated that the expression level of LCK was significantly downregulated as melanoma progressed clinically (p-value< 0.05) ( Figure 5E ). Subsequently, gene set enrichment analyses (GSEA) were carried out in the TCGA cohort and the GEO cohort. The GEO cohort we use here was integrated from 7 separate data, including GSE15605, GSE19234, GSE22155, GSE3189, GSE46517, GSE59455 and GSE65904 ( Supplementary Figure 1 ). The results of the enrichment analysis showed that the gene set in the LCK high expression group was significantly enriched in cancer progression-related pathways and immune component-related pathways ( Figure 5F ). The above results illustrate that LCK plays an important role in remodeling the tumor microenvironment of melanoma.

To explore the specific role played by LCK in the tumor microenvironment of melanoma, we analyzed the immune landscape of LCK in melanoma. To begin, we focused on single-cell sequencing data from melanoma. We integrated single-cell sequencing data from two melanoma cases and further completed the dimensional reduction clustering and annotation ( Supplementary Figures 2A–D , Figure 6A ). We then found that LCK was significantly highly expressed in T cells and could be further used as a marker for T cells ( Figure 6B ). Subsequently, we calculated the ESTIMATE score for each sample in the TCGA cohort through the ESTIMATE algorithm. We further found that ESTIMATE scores were significantly higher in the group with high LCK expression levels than in the group with low LCK expression levels (p-value< 0.05) ( Figure 6C ). Protein-protein interaction network constructed based on LCK showed that LCK plays a key role in the T-cell antigen receptor-linked signal transduction pathway ( Figure 6D ). Immediately after, we classified the samples in the TCGA cohort of melanoma into cold and hot tumors based on the 12 genes extracted from the study of Chunyu Dong et al. ( Supplementary Figure 3 ). We found that LCK expression levels were higher in hot tumors ( Figure 6E ).

Immunotherapy such as immune checkpoint inhibitors have wide application in some solid tumors such as melanoma. However, the issue of patient responsiveness is an obstacle to their effective application. Here, we analyzed the relationship between the expression levels of LCK and the expression levels of the 9 immune checkpoints. The correlation analysis showed that LCK showed a significant positive correlation with immune checkpoints ( Figure 6F ). That is to say, patients with high LCK expression tend to have better immunotherapeutic responses due to high levels of ICPs [31]. Thus, to verify the relationship between LCK and immunotherapy response, we examined the expression levels of LCK in the immunotherapy cohort. Based on 2 immunotherapy cohorts extracted from the TIDE database, we observed that LCK expression levels were higher in the immunotherapy response group ( Figure 6G ). Collectively, in melanoma, LCK is closely linked to the immune components of its microenvironment and is effective in predicting immunotherapeutic response.

We found that LCK, a marker of T cells, plays a crucial role in the immunotherapy of melanoma. Thus, we tried to elucidate the specific mechanism by which LCK acts through cellular interactions. Based on the previously annotated single cell sequencing data, we extracted the melanoma cells and T cells separately. By further dimensional reduction clustering, we obtained 14 subgroups ( Supplementary Figures 4A, B ). Based on the differential expression of genes, we classified melanoma cells into 4 categories ( Supplementary Figure 4C , Figure 7A ). Meanwhile, we classified T cells into Naive CD4⁺ T cells, Naive CD8⁺ T cells, Effector CD8⁺ T cells and Memory CD8⁺ T cells according to the specific expression of CD4, CD8 (CD8A), CD45 (PTPRC), CD197 (CCR7), CD25 (IL2RA) ( Supplementary Figure 4D , Figure 7A ). Subsequently, cell-cell communication was inferred by the ‘CellChat’ package on the basis of the subgroup we annotated. The results showed that there was a strong interaction between melanoma cells and T cells ( Supplementary Figures 4E, F ). Subsequently, we selected LCK-related ligand-receptor pairs left for further analysis. The results suggest that LCK as a ligand acts between the C2 subgroup of melanoma cells and T cell subgroups by binding to the CD8 receptor ( Figures 7B–D ). The binding of LCK to the CD8 receptor drives CD8 close to the TCR and ultimately stabilizes the TCR-MHCp interaction (31). TCR-MHCP interactions promote T cell activation (32). The activation of T cells further increases the response of T cells to pathogens or malignant cells (33). Therefore, we propose that LCK inhibits immune escape of melanoma cells in melanoma by promoting the activation of T cells. Finally, we examined the correlation between the expression levels of LCK and the IC50 of commonly used antitumor drugs. The results showed that the group with high expression of LCK had a lower IC50 ( Figure 7E ). This means that LCK can effectively enhance the inhibitory effect of antitumor drugs.