Landscape of immune-related genes in melanoma
The entirety of this study was accomplished through the following flowchart (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 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 (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 significantly different (p-value < 0.05) (Fig. 1C). Genes with significantly different results 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[26]. The results of mutation profiling of the primary screening genes showed that the mutation rate of all primary screening genes was low (Fig. 1D).
Table 1
A list of core immune-escape related genes.
Gene Symbol
|
ADAR, B2M, BECN1, CALR, ERAP1, FAS, HDAC1, 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
Identification 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 “ClusterProfile” 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 filtered 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 filter 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 finally identified 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 significantly different (Fig. 3C). At last, three hub-independent prognostic factors were obtained by overlapping the independent prognostic factors with significant 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 infiltration 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 significantly 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 infiltration 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 infiltration levels of the six kinds of immune cell in the high- and low-risk groups showed that in the low-risk score group, all six kinds of immune cells had high levels of infiltration (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 infiltration of six kinds of immune cells and tumor purity and results showed that LCK was positively correlated with the infiltration 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 infiltration 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 infiltration 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 infiltration 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 significantly 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 identified LCK as a potential molecular marker for predicting immune checkpoint efficacy. 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 significantly higher in the LCK high-expression group, while the proportion of Macrophages M0 and Macrophages M1 were significantly 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 infiltration levels. While LCK expression was negatively correlated with the infiltration 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 infiltration showed that highe levels of infiltration of T cell CD8 and T cell CD4 memory activated were significantly 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]. From the structural point of view, LCK has the typical organization found in all members of the SRC kinase family comprising an N-terminal site (SRC -homology 4, SH4 domain), a unique region, a SH3 and SH2 domain, a catalytic domain, and a short C-terminal tail[28](Supplementary 6C).