Molecular Proling of Core Immune-Escape Genes Highlights LCK Has Potential to Reshape TME in Melanoma

Background: The game between the immune system of the organism and the tumor is a dynamic course of events. Once the escape phase is reached, the tumor will overcome the limitations of the immune system and the tumor will progress. Objective: This study aimed to determine the potential tumor environment-based prognostic biomarkers related to immunotherapy in melanoma. Methods: 471 tumor samples and 398 normal samples were extracted from The Cancer Genome Atlas and The Genotype-Tissue Expression. Furthermore, a core set of immune-escape genes were collected from previous studies and a set of immune-related genes were obtained from IMMPORT database. Through overlapping these two sets of genes, a set of core immune-escape related genes were identied. In the next place, we conducted a systematic analysis of core immune-escape related genes through The Cancer Genome Atlas cohort and identied independent prognostic factors in melanoma. Through CIBERSORT, ssGSEA algorithm and TIMER database, we explored the potential of independent prognostic factors to reshape the tumor microenvironment. Results: Through Kaplan-Meier as well as univariate cox regression analysis of expression proles and clinical information obtained from the TCGA cohort, we found that high LCK expression was associated with prolonged overall survival. In addition, the expression level of LCK was signicantly correlated with the dysregulation of the inltration level of immune cells in the tumor microenvironment. Conclusions: In this study, we determined LCK as a biomarker that is signicantly associated with tumor environment and has signicant prognostic signicance. Meanwhile, immunotherapeutic approaches targeting LCK in tumor cells may provide a new perspective for the treatment of melanoma. EGFR was associated with shortened overall survival. (C) Correlation analysis of the expression levels of three hub-independent prognostic factors with the inltration levels of six kinds of immune cells and tumor purity. The results showed that the expression level of LCK was highly positively correlated with the inltration level of six kinds of immune cells and highly negatively correlated with tumor purity, revealing the key role of LCK in TME. The pink border represents the absolute value of r correlation coecient is greater than 0.3, and the green border represents the absolute value of r correlation coecient is less than 0.3; the pink border represents the mRNA levels of three hub-independent prognostic factors are positively correlated with the level of six kinds of immune cell inltration and tumor purity, and the green border represents the mRNA levels of three hub-independent prognostic factors are positively correlated with the level of six kinds of immune cell inltration and tumor purity. negative correlation. (D) Correlation analysis of the expression levels of three hub-independent prognostic factors with the inltration levels of 29 kinds of immune cells. Except for mast cells, LCK expression was highly positively correlated with the level of inltration of the remaining 28 kinds of immune cell types. (C) (D) The difference of proportion of between and represents a of cells; represents low T cell CD8 immune cell fraction was in the group, while Macrophages M0 and M2 (E) Correlation analysis of

The immune system has a crucial role in the development and treatment of cancer. Adaptive immunity can prevent or constrain cancer through immunosurveillance, whereas innate immunity and in ammation tend to promote tumorigenesis and malignant development of nascent cancer [6]. Cancer immunoediting is the process by which the immune system both constrains and promotes tumor development, and it proceeds through three phases: elimination, equilibrium and escape [7]. Cancer is a highly heterogeneous disease involving complex interactions between the malignant cells and the tumor microenvironment (TME) [8]. Cancer cells have evolved multiple mechanisms to evade immune recognition or to modulate immune cell function [9]. TME, a mixture that consists of mesenchymal cells, tumor-in ltrating immune cells (TIICs), endothelial cells, extracellular matrix molecules and in ammatory mediators [10]. Increasing evidence demonstrated the importance of the TME in the tumor development [11]. Collaborative interactions between cancer cells and their supporting cells contributed to the malignant phenotypes of cancer, such as immortal proliferation, resisting apoptosis, and evading immune surveillance. Therefore, the TME signi cantly in uences therapeutic response and clinical outcome in cancer patients [12,13]. TME is closely correlated to tumor initiation, progression and prognosis [14]. The composition of the tumor-in ltrating population re ects the mechanisms underlying anticancer immune responses and can help identify novel prognostic signatures [8]. thereby assessing the immune in ltration of tumor microenvironment species [15]. Furthermore, single cell gene set enrichment analysis (ssGSEA) algorithm is a new algorithm developed on the basis of GSEA, rst proposed in 2009 by David A. Barbie et al. This was achieved through a single-sample extension of GSEA, designed to allow researchers to de ne an enrichment score that would allow one to obtain the level of in ltration of 29 immune cell species in each sample [16]. In this study, we obtained a set of genes related to immune escape from an article and identi ed differentially expressed proteins regulated by immune-escape related genes through the cbiportal website, identi ed protein-coding genes highly associated with prognosis by constructing a prognostic model, and nally described the role played by independent prognostic factors in remodeling TME using CIBERSORT and ssGSEA.

Data collection
Gene expression pro les, clinical phenotype information and somatic mutation data of melanoma were extracted from The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov) database. Melanoma sample information was obtained from TCGA and Genotype-Tissue Expression Project (GTEx). The expression pro le data are presented in TPM form and are normalized by log2(x + 0.001). TCGA data were downloaded through the "TCGAbiolinks" package through R project and the above pre-processing of gene expression pro les was carried out in R project.
Identi cation and characterization of immune-related genes in melanoma A core set of immune-escape related gene were extracted from papers [17]and interference with these genes alone will increase the sensitivity or resistance of cancer cells to Cytotoxic T lymphocyte (CTL)mediated toxicity. Furthermore, a set of immune-related genes was acquired from the IMMPORT (https://www.immport.org/) database which were curated through functions and gene ontology terms. To be next, overlapping these two sets of gene to obtain core immune-escape related genes. In addition, a combination of univariate cox regression analysis and Kaplan-Meier survival analysis was utilized to evaluate the prognostic value of the intersecting genes, while differential expression analysis of crossover genes in normal and tumor tissues was conducted in order to determine primary screening genes. Subsequently, we visualized the somatic mutation data in Mutation Annotation Format (MAF) using the "maftools" R package, which provides a large amount of commonly used analysis and visualization modules in cancer genomic studies [18].
Identi cation of differentially expressed proteins and their enrichment pathways, analysis of proteinprotein interactions.
The cBioPortal for Cancer Genomics(http://cbioportal.org) provides a Web resource for exploring, visualizing, and analyzing multidimensional cancer genomics data [19]. The protein expression changes (both unphosphorylated and phosphorylated proteins) in melanoma were analyzed in cBioPortal through Reverse Phase Protein Array (RPPA) [20]. Proteins with signi cant differences (-log10 P-value > 1.3) are indicated by dots of different colors (The under-expression group is indicated in green and the overexpression group is indicated in red). The R package "ggplot2" and "ClusterPro ler" were then used to perform enrichment analysis of these protein-coding genes and visualize the results, which contains Gene Ontology (GO) enrichment analysis of Cellular component (CC), Biological process (BP) and Molecular function (MF) as well as Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis. STRING database (https://string-db.org/cgi/input.pl) can collect and integrate known and predicted protein-protein association data [21]. In this study, we used STRING database to construct Protein-Protein Interaction (PPI) network with minimum required interaction score 0.4 and disconnected nodes in the network were hide. Also, we exhibited the nodes with the top 60 highest connectivity in the network. At last, network was visualized through Cytoscape software and three important subnetworks were selected by MCODE which is a plugin within Cytoscape with a degree cut-off = 2, node score cut-off = 0.2, k-core = 2, and max Depth = 100 as the criterion.

Risk modeling and screening of independent prognostic factors
The melanoma patients from TCGA 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. The R package "survival" and "UpSetR" was used to perform univariate cox regression analysis to screen for protein-coding genes for signi cant association with overall survival in the TCGA-SKCM cohort. Protein-coding genes with p-value less than 0.05 were used for further analysis, and the t was reduced by least absolute shrinkage and selection operator (LASSO) regression analysis to obtain fewer representative eigenvalues. 1000-round cross-validation for tuning parameter selection was used to prevent over tting and the partial likelihood deviance met the minimum criteria [22]. The screened protein-coding genes were then subjected to multivariate cox proportionalhazard analysis through R package "survival" in R project to identify independent prognostic factors (pvalue less than 0.05). In addition, three hub-independent prognostic factors were determined by overlaying the nodes with the top 60 connectivity in PPI with independent prognostic factors.
Validation and analysis of risk prognostic models Risk score were calculated for each patient sample through the coef values from the multivariate cox proportional-hazards analysis results and the expression levels of independent prognostic factors in the results of multivariate analysis with the following formula: Patients were then divided into high-and low-risk groups according to the median risk score, and a risk score model was developed. Then, ROC curves are created through the R package "pROC" to verify the accuracy of the risk model, with larger Area Under Curve (AUC) representing more accurate models. This risk model is then subjected to a series of analyses based on the established risk model. First of all, risk factor analysis was conducted through the online mapping site HIPLOT (https://hiplot/com.cn) to evaluate the association between patient survival status and expression level of independent prognostic factors and risk scores. Furthermore, we divided the sample into high-and low-risk score groups and assessed the ability of the risk scores to predict patient prognosis by Kaplan-Meier survival analysis. At last, through TIMER, a web server for comprehensive analysis of tumor-in ltrating immune cells, data on the in ltration level of six kinds of immune cells were downloaded [23]. Differential in ltration level of six kinds of immune cell in high-and low-risk score subgroups was analyzed.
The expression, survival and immune microenvironment analysis of Hub-independent prognostic factors Prognostic value of all hub-independent prognostic factors and differential expression in tumor and normal tissues were analyzed in the TCGA cohort. Moreover, to further demonstrate that the accuracy of the above results, Gene Expression Pro ling Interactive Analysis (GEPIA, http://gepia.cancer-pku.cn/), a web server for cancer and normal gene expression pro ling and interactive analysis [24], was utilized to validated above results. Furthermore, the association between 6 kinds of immune cells and hubindependent prognostic factors was analyzed through TIMER database. Similarly, we quanti ed the level of in ltration of 29 kinds of immune cell types through the ssGSEA algorithm in the R package "GSVA" and then assessed the relationship between the in ltration level of 29 kinds of immune cell and the mRNA expression of three hub-independent prognostic factors.
The relationship between LCK with the proportion of (tumor-in ltrating immune cells) TICs In this study, CIBERSORT algorithm was used to evaluate the proportion of 22 kinds of immune cells in each sample, which allowed us to assess the correlation between the mRNA expression of LCK and the level of in ltration of immune components. Only samples with p-value less than 0.05 were included in the next further analysis.

Statistical analysis
Analysis of variance between two groups was performed by Student's t-test, and analysis of variance between multiple groups (n = > 3) was performed by one-way ANOVA. Prognostic survival analysis was conducted through Kaplan-Meier as well as univariate and multivariate risk proportional regression models. p-values were calculated through log-rank test. p-values less than 0.05 were considered to be signi cantly different. Correlation coe cients were calculated using Spearman, and correlation coe cients with absolute values greater than 0.3 were considered to be correlated. The above statistical analyses were performed in R project as well as in PRISM 7.0 software.

Landscape of immune-related genes in melanoma
The entirety of this study was accomplished through the following owchart (Supplementary 1). A set of immune-escape related genes and a set of immune-related genes were extracted from the literature 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 1) were determined as subjects for further analysis (Fig. 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. Results of Kaplan-Meier survival analysis of core immune-escape related genes as well as univariate cox analysis showed that except for CALR, TNFRSF1A, HDAC1, JAK1 and TFRC, which were not signi cantly different, and MAPK1, which was an unfavorable prognostic factor, the rest of all the core immuneescape genes were favorable prognostic factors, and all of them were signi cantly different (Fig. 1B,  Supplementary 2). Furthermore, the results of differential expression analysis of core immune-escape related genes in tumor and normal tissues showed that 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 signi cantly different (p-value < 0.05) (Fig. 1C). Genes with signi cantly different results in both expression difference analysis and prognostic analysis were initially screened out for inclusion in the next analysis. Somatic mutation pro les of 467 melanoma patients downloaded from TCGA database were analyzed and visualized via the "maftools" R package [26]. The results of mutation pro ling of the primary screening genes showed that the mutation rate of all primary screening genes was low (Fig. 1D).  IFNAR1, IFNAR2, IRF1, IRF9, IFNGR1, IFNGR2, IKBKG,  IKBKB, JAK1, JAK2, MAPK1, PDLA3, PSMB8, SOCS1, STAT1, TAP1, TAP2, TAPBP, TBK1, TFRC,  TGFBR2, TNFAIP3, TNFRSF1A, TNFRSF1B   Table 2. Representative Gene Ontology results of these three subnetworks.
The differentially expressed proteins under the regulation of primary screening genes are mainly enriched in multiple immune-related pathways and functional Identi cation of differentially expressed proteins regulated by primary screening genes through the RPPA module of the cBioPortal website and 164 differentially expressed proteins were obtained in total. Subsequently, GO and KEGG enrichment analysis of these protein-coding genes was conducted through the "ClusterPro le" package in R project. The pathway enrichment results showed that these protein-coding genes were mainly enriched in EGFR tyrosine kinase inhibitor resistance, ErbB signaling pathway, PI3K-Akt signaling pathway and some immune-related pathways (including T cell receptor signaling pathway, Natural killer cell mediated cytotoxicity and so on) ( Fig. 2A). GO enrichment analysis showed that these protein-coding genes were mainly enriched in the regulation of apoptotic signaling pathway and regulation of mitochondrion organization in BP, chromosome, telomeric region and mitochondrial outer membrane in CC, phosphatase binding and ubiquitin-like protein ligase binding in MF (Fig. 2B). Through STRING database, protein-protein interaction network was constructed, which contains 188 nodes and 3670 edges in summary. After that, we ltered out the top 60 nodes in the entire network in terms of connectivity (Fig. 2C). Afterwards, we used the MCODE plugin in Cytoscape software in order to lter out the three important sub-networks of the network (Fig. 2D). Finally, GO enrichment analysis was performed on the nodes in each of the three subnetworks respectively, and the results showed that the main enrichment was in terms related to immunization ( Table 2).
Construction and evaluation of the melanoma-related risk signature which are predictive of prognosis The samples are divided into a TRS and a TES according to an approximate ratio of 1:1 (The TRS has 228 samples and the TES has 226 samples). The result of univariate cox regression analysis of TRS showed that a total of 15 prognostic factors were determined, which means changes in the expression levels of these 15 protein-coding genes affect the prognosis of patients. In addition that, except CDKN1B, LCK, and RICTOR which were favorable prognostic factors, the rest of these prognostic factors were associated with reduced overall survival (the p-value of the above results were less than 0.05) (Fig. 3A). Subsequently, prognostic factors were subjected to LASSO regression analysis and 10 representative protein-coding genes were determined (Fig. 3B). Multivariate cox regression analysis was conducted within these 10 representative protein-coding genes obtained from LASSO regression analysis and six independent prognostic factors were nally identi ed which were related to prognosis in melanoma. In summary, all four independent prognostic factors were associated with reduced overall survival, except for LCK and RICTOR, which were favorable factors, although the results for RICTOR were not signi cantly different ( Fig. 3C). At last, three hub-independent prognostic factors were obtained by overlapping the independent prognostic factors with signi cant differences in the multivariate results and the top 60 nodes with the highest connectivity in the protein interactions network (Fig. 3D).
Risk-prognosis models constructed by independent prognostic factors may prolong overall survival by affecting immune cell in ltration in TME Through analyzing the independent prognostic factors in the TRS, a risk prognostic model was constructed in the TRS and TES, using which the prognosis of patients could be effectively predicted. In the TRS and TES, ROC curves were established to validate the accuracy of the risk model, with AUC values equal to 0.683 and 0.657, respectively, indicating that the model has good accuracy (Fig. 4A). After that, the results of the risk factor analysis showed that the number of patient deaths clustered signi cantly as the risk score increased, indicating that the risk score can effectively predict the prognosis of patients. In both TRS and TES, LCK and RICTOR were higher expressed in the low-risk group, while FOXM1, KIT, EGFR, SMAD1 were higher expressed in the high-risk group (Fig. 4B). Survival analysis in the high-and low-risk score groups showed that the high-risk group was associated with lower overall survival (p-value is less than 0.05) (Fig. 4C). The in ltration levels of six kinds of immune cells were extracted from the TIMER database, including B cell, CD4 T cell, CD8 T cell, Neutrophil, Macrophage and Dendritic. Analysis of the differences in the in ltration levels of the six kinds of immune cell in the highand low-risk groups showed that in the low-risk score group, all six kinds of immune cells had high levels of in ltration (p-value is less than 0.05) (Fig. 4D).
Aberrant expression and prognostic value analysis of hub-independent Prognostic factors in melanoma and LCK might serve as a promising indicator for remodeling TME The results of differential expression of three hub-independent prognostic factors in normal and tumor tissues through the TCGA cohort showed that LCK was highly expressed in tumor tissues, while KIT and EGFR were highly expressed in normal tissues (Fig. 5A). In addition, the results of Kaplan-Meier survival analysis of the three hub-independent prognostic factors, also derived from TCGA cohort, showed that high-expression of LCK was associated with prolonged overall survival, while high expression of KIT and EGFR was associated with reduced overall survival (Fig. 5B). The results of the above analysis were reproduced through the GEPIA database (Supplementary 4). Furthermore, we analyzed the correlation between the expression of the three hub-independent factors and clinical parameters, including lymph nodes, tumor topography, pathologic stage, event, gender and metastasis (Supplementary 5). The results of analysis showed that LCK expression was negatively correlated with the size of the primary tumor (tumor topography) (p-value = 3.6e-05) and higher degree of metastasis correlates with lower LCK expression (p-value = 0.0024). LCK expression was higher in surviving patients, while KIT and EGFR expression was higher in deceased patients. Through TIMER database, we analyzed the relationship of the mRNA expression of three hub-independent prognostic factors with the level of in ltration of six kinds of immune cells and tumor purity and results showed that LCK was positively correlated with the in ltration level of the six kinds of immune cells and negatively correlated with tumor purity (p-value is less than 0.05) (Fig. 5C). Besides, correlation analysis of the in ltration levels of 29 immune cell calculated by the ssGSEA algorithm with the expression levels of three hub-independent prognostic factors showed that LCK is highly positively correlated with the in ltration level of 28 kinds of immune cells other than mast cells (Fig. 5D), indicating that LCK play an important role in melanoma occurrence and progression partly because of altering the level of in ltration of immune cells in TME.

Correlation analysis of the LCK to common ICPs and to the proportion of TICs
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. Therefore, in this study, we analyzed the relationship between the three hub-independent prognostic factors and 12 kinds of ICPs. Correlation analysis showed that LCK was signi cantly and positively correlated with 10 kinds of ICPs other than IL6 and CD276 (Fig. 6A). Furthermore, high expression of ICPs was observed in high LCK expression group (Fig. 6B). The results demonstrated that patients with high LCK expression tended to have a better immunotherapy response because of the high levels of ICPs [25]. In light of the above series of analyses, we identi ed LCK as a potential molecular marker for predicting immune checkpoint e cacy. After that, we divided the samples into two phenotypic cohorts according to high-and low-LCK expression and used the c2.cp.kegg.v7.4.symbols.gmt gene set for GSEA analysis, showing that the LCK high expression phenotype was mainly enriched in antigen-processing and presentation, cell adhesion molecules cams, chemokine signaling pathway, cytokine-cytokine receptor interaction, FC-GAMMA-R mediated phagocytosis, leukocyte transendothelial migration and natural killer cell mediated cytotoxicity. While LCK low expression phenotype was mainly enriched in ribosome (Fig.  6C). The results of CIBERSORT analysis showed that T cells CD8, macrophages M0, macrophages M1 and macrophages M2 were relatively high in all samples within 22 kinds of immune cell. ( Supplementary  6A). Furthermore, analysis of the relationship between LCK expression and 22 kinds of immune components showed that the proportion of T cell CD8 was signi cantly higher in the LCK high-expression group, while the proportion of Macrophages M0 and Macrophages M1 were signi cantly lower in the LCK low-expression group (Fig. 6D). Moreover, the results of the correlation analysis between immune components and LCK expression showed that LCK expression was positively correlated with T cell regulatory (Tregs), T cell CD4 memory activation, and T cell CD8 in ltration levels. While LCK expression was negatively correlated with the in ltration level of Macrophages M0, Macrophages M2, NK cells resting (Fig. 6E).

The Kaplan-Meier survival analysis based on the level of T cell CD8 and T cell CD4 memory activated immune cell in ltration showed that highe levels of in ltration of T cell CD8 and T cell
CD4 memory activated were signi cantly associated with increased overall survival (Supplementary 6B).
The above results suggested that LCK was associated with antitumor immunity in melanoma, which partially explained the association of LCK with better prognosis [27].

Discussion
Cancer cells must acquire phenotypic changes that allow them to evade recognition and destruction by effector cells of the immune system such as CTLs [17]. Factors that likely relate to immune escape include lack of strong cancer antigens or epitopes recognized by T cells, minimal activation of cancerspeci c T cells, poor in ltration of T cells into tumors, downregulation of major histocompatibility complex on cancer cells, and immunosuppressive factors and cells in the tumor microenvironment [29]. We identi ed a set of immune-escape related genes from the literature and overlapped them with immune-related genes to obtain core immune-escape genes. Through a series of comprehensive analysis of these immune-escape related genes, the results suggest that LCK could be a new potential target that could predict the overall survival time of melanoma patients. [25] and its high expression was signi cantly associated with elevated ICPs. KEGG enrichment analysis of 164 protein-coding genes, including LCK, showed that LCK is enriched in multiple immune-related pathways, including PD-L1 expression and the PD-1 checkpoint pathway in cancer, T-cell receptor signaling pathway, natural killer cell-mediated cytotoxicity, Th17 cell differentiation and Th1 and Th2 cell differentiation. GO enrichment analysis showed that LCK is involved in T cell activation, SH2 domain binding, Protein tyrosine kinase activity, Regulation of protein kinase B signaling, Membrane region, Phosphatase binding, Protein Cterminal binding, Protein phosphorylated amino acid binding, and Phosphotyrosine residue binding. In addition, the results of the CIBERSORT analysis, which calculates the proportion of immune components in the samples, showed a higher proportion of antitumor immune cells (T cell CD8) and a lower proportion of protumor immune cells (Macrophages M0 and Macrophages M2) in the LCK high expression group. T cell CD8 is an essential component of tumor-in ltrating lymphocytes (TIL) and TME [30]. The roles and prognostic value of TILs have been extensively studied in various malignant tumors. TILs are known to have critical effects on the survival of malignant melanoma and act as prognostic markers [31][32][33][34]. TIMER website analysis showed that LCK expression was positively regulate the quantity and differentiation of immune cells [35]. Immune cells are the main components of TME, and their number and status play a critical role in the progression of tumor development, invasion and metastasis [36][37][38]. Therefore, elucidating the in ltrating immune cells of TME may help to elucidate the potential molecular mechanisms of LCK involvement in melanoma [27]. Indeed, the uctuating expression levels of LCK do affect the level of in ltration of immune cells in TME, which means the high level of LCK expression measured by bulk-seq is indeed a key factor that can reshape TME.
LCK (lymphocyte-speci c protein tyrosine kinase) belongs to the SRC family of tyrosine kinases and has been best studied in the context of T-cell function and signaling as well as lymphocytic leukemia of the Bcell lineage [39]. LCK is mainly expressed in T cells, NK cells, brain, B cells [39,40] and plays a vital role in various cellular processes such as cell cycle control, cell adhesion, motility, proliferation and differentiation [41]. This encoded protein is a key signaling molecule for selection and maturation during T cell development. LCK was strongly expressed in lymphocytes. However, tumor cells also expressed LCK, albeit to a lesser extent [42]. In recent years, it is reported that LCK has emerged as one of the key molecules regulating T cell function, and studies using knockout LCK mice or LCK-de cient T cell lines surface that LCK regulates signal initiation, T cell development and T cell homeostasis and also can enhan or inhibit BCR signaling. Patients with LCK de ciency frequently present with immune dysregulation and autoimmunity. Over expression of LCK leads to large number of other diseases like cancer, asthma, diabetes 1, rheumatoid arthritis, psoriasis, systemic lupus erythematosus, in ammatory bowel diseases (crohn's disease and ulcerative colitis), organ graft rejection, atherosclerosis, hypersensitivity reactions, polyarthritis, dermatomyositis [43,44]. Since LCK is involved in T cell proliferation and differentiation, therefore, new small molecules with LCK inhibitory activity can be of great relevance to treat T cell mediated diseases. It has been reported that increased LCK expression sometimes leads to colorectal cancer [45]. In addition, another report stated that inhibition or downregulation of LCK led to apoptosis in CLL cell lines [46]. Therefore, the application of LCK inhibitors could be an important strategy for the treatment of certain cancers [47]. However, it has also been reported that high expression of LCK is associated with increased cumulative survival in melanoma patients [48].
In another cohort survival analysis report on LCK-de cient CLL disease models in mice, results showed that LCK-KO control mice exhibited a signi cantly shorter median survival time compared to wild-type healthy mice over an observation period of 350 days due to T-cell de ciency and the resulting immunode ciency [49]. Likewise, in our present study, elevated LCK expression levels in the TCGA cohort had an independent prognostic value in melanoma, and results demonstrating that upregulation of LCK expression plays a defensive role in the developmental progression of melanoma. This was reproduced in the GEPIA and the human atlas databases (HPA). LCK functions primarily in lymphocytes and is involved in transduction from the T-cell receptor complex to the nucleus, and this speci c expression and function in immune cells may partially explain the phenomenon that high LCK expression is often associated with extended overall survival. On the other hand, RNA-seq is actually a kind of bulk-seq, in which various cells are mixed, so the highly expressed genes measured by RNA-seq in tumor tissues are likely to be molecules highly expressed on immune cells. In addition, since LCK plays a role in cancer cell signaling as well as in T-cell function, it will be necessary to de ne therapeutic strategies to selectively target LCK in tumor cells without impairing the responses of tumor in ltration lymphocytes. This is a critical issue common to other kinase inhibitors targeting signaling molecules expressed in both cancer and immune cells (e.g., B-Raf, AKT, mTOR inhibitors) [39].

Conclusion
In the present study, we conclusively identi ed LCK as a molecule that can remodel TME and is associated with prolonged overall survival in melanoma through a systematic analysis of immuneescape related genes, which suggest that LCK is a potential molecular marker which can predict the overall survival in melanoma. In addition, LCK targeting tumor cells will also be a promising therapeutic approach.    accurate. (B) Risk factor analysis in TRS and TES, the samples are divided into two groups of high-and low-risk according to the risk score, with green dots representing surviving samples and red dots representing dead samples. As the risk score increases, the number of dead samples also increases. The heatmap shows the expression of the six factors in the multivariate results in the high-and low-risk groups. Green area represents low expression; red area represents high expression. As the risk score increased, the mortality sample also increased, indicating that the risk score model does predict patient prognosis. (C) Survival curves based on high-and low-risk score for the training and test sets. In both the training and test sets, low-risk scores were associated with elevated overall survival. (D) Differential expression analysis of six immune cells in high-and low-risk groups, red boxes represent the high-risk group, blue boxes represent the low-risk group. Low-risk subgroups have higher in ltration levels of immune cells in ltrating into the TME.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.