Figure 1.
Workflow of this study.
The tumor microenvironment plays a crucial role in melanoma. In this study, the abundance of immune cells in melanoma samples was assessed and analyzed using single sample gene set enrichment analysis (ssGSEA), and the predictive value of immune cells was assessed using univariate COX regression analysis. The Least Absolute Shrinkage and Selection Operator (LASSO)-Cox regression analysis was applied to construct an immune cell risk score (ICRS) model with a high predictive value for identifying the immune profile of melanoma patients. The pathway enrichment between the different ICRS groups was also elucidated. Next, five hub genes for diagnosing the prognosis of melanoma were screened by two machine learning algorithms, LASSO and random forest. The distribution of hub genes in immune cells was analyzed on account of Single-cell RNA sequencing (scRNA-seq), and the interaction between genes and immune cells was elucidated by cellular communication. Ultimately, the ICRS model on account of two types of immune cells (Activated CD8 T cell and Immature B cell) was constructed and validated, which can determine melanoma prognosis. In addition, five hub genes were identified as potential therapeutic targets affecting the prognosis of melanoma patients.
Citation: Linqian Guo, Qingrong Meng, Wenqi Lin, Kaiyuan Weng. Identification of immune subtypes of melanoma based on single-cell and bulk RNA sequencing data[J]. Mathematical Biosciences and Engineering, 2023, 20(2): 2920-2936. doi: 10.3934/mbe.2023138
| [1] | Jianpei Hu, Zengnan Mo . Dissection of tumor antigens and immune landscape in clear cell renal cell carcinoma: Preconditions for development and precision medicine of mRNA vaccine. Mathematical Biosciences and Engineering, 2023, 20(2): 2157-2182. doi: 10.3934/mbe.2023100 |
| [2] | Chenyang Jiang, Weidong Jiang . AGTR1, PLTP, and SCG2 associated with immune genes and immune cell infiltration in calcific aortic valve stenosis: analysis from integrated bioinformatics and machine learning. Mathematical Biosciences and Engineering, 2022, 19(4): 3787-3802. doi: 10.3934/mbe.2022174 |
| [3] | Xiaowei Zhang, Jiayu Tan, Xinyu Zhang, Kritika Pandey, Yuqing Zhong, Guitao Wu, Kejun He . Aggrephagy-related gene signature correlates with survival and tumor-associated macrophages in glioma: Insights from single-cell and bulk RNA sequencing. Mathematical Biosciences and Engineering, 2024, 21(2): 2407-2431. doi: 10.3934/mbe.2024106 |
| [4] | Kaiyu Shen, Shuaiyi Ke, Binyu Chen, Tiantian Zhang, Hongtai Wang, Jianhui Lv, Wencang Gao . Identification and validation of biomarkers for epithelial-mesenchymal transition-related cells to estimate the prognosis and immune microenvironment in primary gastric cancer by the integrated analysis of single-cell and bulk RNA sequencing data. Mathematical Biosciences and Engineering, 2023, 20(8): 13798-13823. doi: 10.3934/mbe.2023614 |
| [5] | Yi Shi, Xiaoqian Huang, Zhaolan Du, Jianjun Tan . Analysis of single-cell RNA-sequencing data identifies a hypoxic tumor subpopulation associated with poor prognosis in triple-negative breast cancer. Mathematical Biosciences and Engineering, 2022, 19(6): 5793-5812. doi: 10.3934/mbe.2022271 |
| [6] | Yu Jiang, Lijuan Lin, Huiming Lv, He Zhang, Lili Jiang, Fenfen Ma, Qiuyue Wang, Xue Ma, Shengjin Yu . Immune cell infiltration and immunotherapy in hepatocellular carcinoma. Mathematical Biosciences and Engineering, 2022, 19(7): 7178-7200. doi: 10.3934/mbe.2022339 |
| [7] | Jun Wang, Mingzhi Gong, Zhenggang Xiong, Yangyang Zhao, Deguo Xing . Immune-related prognostic genes signatures in the tumor microenvironment of sarcoma. Mathematical Biosciences and Engineering, 2021, 18(3): 2243-2257. doi: 10.3934/mbe.2021113 |
| [8] | Rui-zhe Zheng, Jin Xing, Qiong Huang, Xi-tao Yang, Chang-yi Zhao, Xin-yuan Li . Integration of single-cell and bulk RNA sequencing data reveals key cell types and regulators in traumatic brain injury. Mathematical Biosciences and Engineering, 2021, 18(2): 1201-1214. doi: 10.3934/mbe.2021065 |
| [9] | Jiaxin Luo, Lin Wu, Dinghui Liu, Zhaojun Xiong, Linli Wang, Xiaoxian Qian, Xiaoqiang Sun . Gene regulatory network analysis identifies key genes and regulatory mechanisms involved in acute myocardial infarction using bulk and single cell RNA-seq data. Mathematical Biosciences and Engineering, 2021, 18(6): 7774-7789. doi: 10.3934/mbe.2021386 |
| [10] | Rongxing Qin, Lijuan Huang, Wei Xu, Qingchun Qin, Xiaojun Liang, Xinyu Lai, Xiaoying Huang, Minshan Xie, Li Chen . Identification of disulfidptosis-related genes and analysis of immune infiltration characteristics in ischemic strokes. Mathematical Biosciences and Engineering, 2023, 20(10): 18939-18959. doi: 10.3934/mbe.2023838 |
The tumor microenvironment plays a crucial role in melanoma. In this study, the abundance of immune cells in melanoma samples was assessed and analyzed using single sample gene set enrichment analysis (ssGSEA), and the predictive value of immune cells was assessed using univariate COX regression analysis. The Least Absolute Shrinkage and Selection Operator (LASSO)-Cox regression analysis was applied to construct an immune cell risk score (ICRS) model with a high predictive value for identifying the immune profile of melanoma patients. The pathway enrichment between the different ICRS groups was also elucidated. Next, five hub genes for diagnosing the prognosis of melanoma were screened by two machine learning algorithms, LASSO and random forest. The distribution of hub genes in immune cells was analyzed on account of Single-cell RNA sequencing (scRNA-seq), and the interaction between genes and immune cells was elucidated by cellular communication. Ultimately, the ICRS model on account of two types of immune cells (Activated CD8 T cell and Immature B cell) was constructed and validated, which can determine melanoma prognosis. In addition, five hub genes were identified as potential therapeutic targets affecting the prognosis of melanoma patients.
Skin Cutaneous Melanoma (SKCM) is a deadly and invasive kind of skin cancer that develops when melanocytes in the skin epidermis' basal layer undergo malignant transformation. The incidence of melanoma is rising faster in developed countries. There will be approximately 99,780 new patients of melanoma in the United States in 2022, nearly double the 53,600 cases in 2002 [1]. Environmental, immunological, and genetic factors contribute to melanoma development. The primary treatment for non-metastatic melanoma is surgery. However, not all patients with the disease are candidates for it due to anatomical location, number of lesions, and recurrence rate following surgical resection. For the aforementioned circumstance, the response rate has been enhanced by the development of new medications with specific targeting and immunotherapy. Effective melanoma treatment options include immunotherapy and targeted therapy. Currently, immune checkpoint medications that target CTLA-4 and PD-1 monoclonal antibodies are routinely utilized throughout the world to treat melanoma [2]. Relevant clinical studies have demonstrated good patient response to immunotherapy, such as immune checkpoint inhibitors. Unfortunately, many patients with the disease develop immune resistance to immunotherapy, which may be due to the complicated interactions between tumor cells and the tumor microenvironment [3]. The tumor microenvironment is heterogeneous and consists of various types of cells and extracellular matrix with unique Spatio-temporal interactions between these cellular components that influence the behavior of the tumor [4]. During melanoma development, the proliferation and apoptosis of tumor cells are determined by the activity of immune cells [5]. In the background of cancer immunotherapy, intratumoral heterogeneity influences the therapeutic effectiveness of anticancer medications. The tumor microenvironment (TME), which is heterogeneous and has unique spatiotemporal interactions between all cellular constituents, can be rationally explained as the source of therapeutic pressure and cancer adaptability [6]. It has been shown that tumor immune subtypes can be determined on account of the combined abundance of immune cells, which will help select suitable therapeutic strategies [7]. Therefore, constructing an immune cell risk score that can predict prediction is crucial for melanoma prevention and treatment.
ScRNA-seq technology has recently made advancements that have made it easier to comprehend the heterogeneity of tumor cells. ScRNA-seq analysis is unquestionably a valuable method for learning the characteristics of various cell types within cells and is a critical component of single-cell sequencing technology [8]. For instance, a recent study found that scRNA-seq can analyze the genetic signature of drug-resistant clusters, find clusters that cause drug resistance, and assess cell clusters in tumor tissue [9]. ScRNA-seq can enhance the interpretation of how vital cellular constituents co-stimulate tumor emergence behavior and describe tumors' cellular composition [10]. This investigation field is also an essential part of cancer investigation. However, applying scRNA-seq in this field needs further in-depth studies.
The RNA-Seq TPM (Transcripts Per Kilobase of exon model per Million mapped reads) data and clinical data of SKCM samples were downloaded from the TCGA database (https://portal.gdc.cancer.gov/https://www.gtexportal.org/home/index.html). The corresponding information for melanoma was downloaded from the Gene Expression Omnibus (GEO) [11] (https://www.ncbi.nlm.nih.gov/) database for scRNA-seq data: GSE72056 and bulk RNA-seq data: GSE65904. The dataset of TCGA was utilized as the training set, and the GSE65904 cohort was utilized as the verification set.
The ssGSEA method was employed to calculate the abundance of 23 immune cells in skin cutaneous melanoma according to the expression data of tumor samples [12]. Prognosis-relevant immune cells were screened on account of univariate COX regression analysis, followed by constructing the ICRS with predictive value on account of prognosis-relevant immune cells using the Lasso [13]. The immune cell risk score was calculated as follows:
| ICRS=n∑i=1Coef_i×Xi |
where Coef represents the coefficient of prognostic immune cells and X denotes the abundance of predictive immune cells. The "glmnet" package was utilized to perform LASSO-COX regression analysis to remove excess prognostic immune cells in order to build an immune cell risk score subsequently. The entire sample was classified as different ICRS on account of the median risk score of the training cohort. Kaplan-Meier survival curves were plotted to determine the prognosis value of the immune cell risk score for melanoma. The "ESTIMATE" algorithm was utilized to assess the immune score, stromal score, and ESTIMATE score for each patient in the melanoma cohort [14]. CIBERSORT is a generalized deconvolution algorithm for quantifying the cellular fraction of a large number of tissues on account of gene expression matrices to evaluate the proportion of immune cells of 22 phenotypes. To confirm the precision of the inverse convolution algorithm, the threshold was set to a P value < 0.05.
GSVA switched the analysis object from genes to gene sets for pathway-level differential analysis. GSVA portrays the overall differences between different immune cell risk scores on account of the training cohort. The Wilcoxon rank sum test was next utilized to calculate differential genes between different ICRS. |log2FoldChange| > 1 and P-Value < 0.05 were utilized as significance criteria to determine the differential genes. Gene Ontology (GO) enrichment analysis and Genomes (KEGG) pathway analysis of differential genes were performed using the "clusterProfiler" package of R software.
Two machine learning methods were employed to bolster hub genes affecting melanoma prognosis. The LASSO is utilized to build a penalty function to obtain a more refined model by using the "glmnet" package. It is a biased estimation for data with multiple covariances that improve statistical models' prediction accuracy and understandability. The package "randomForest" [15] was utilized to build a random forest model to filter DEGs. In the random forest classifier, the number of random seeds and decision trees is set to 123,456 and 300, respectively, to obtain a stable model error and high accuracy. The differential genes with utmost significance were screened as essential genes affecting melanoma prognosis for subsequent analysis. For the accuracy of the results, the overlapping genes of LASSO and RF are regarded as hub genes affecting melanoma's prognosis [16].
The scRNA-seq was downloaded from the GEO database, and the study included 4645 cells from 19 patients. Cells with more than 10,000 genes and less than 200 genes were discarded. These cells did not include mitochondrial genes. Ultimately, a total of 4612 cells were included in this study for follow-up analysis. First, we normalized the merged data by log-normalization and identified the first 2000 highly variable genes by the "FindVariableFeatures" function. A principal component analysis (PCA) was subsequently performed to select significant main components (PCs) using the "JackStraw" and "ScoreJackStraw" functions. The "FindAllMarkers" function is utilized to recognize characteristic genes for each cell population [17]. The "RunTSNE" function is then utilized for cell clustering and visualization analysis of TSNE. T-SNE is an algorithm that converts multidimensional data into two-dimensional or three-dimensional data by using accountable flow theory and topological algorithm while maintaining the original data distribution and similarity to achieve the effect of dimensionality reduction. Marker genes are subsequently annotated using the "singleR" package and corrected according to their features using CellMarker.
Cellular communication between cell types is assessed by the "CellChat" package [18]. Specifically, Cellular communication is the process by which a portion of cells in vivo sends signals and the target cells receive them and translate them into changes in cellular function. The mechanism of communication between individual cell types was obtained by comparing the differences in ligand and receptor gene expression between cells of different sample groups across cell types. In this process, the ligand is considered the output signal and the receptor is considered the input signal. The ligand-receptor interaction score reflects the possible interactions between cells.
We downloaded 812 normal skin samples from the GTEx database for joint analysis, with 463 melanoma samples from the TCGA database. Boxplot visualizes the differential expression levels of hub genes between tumor and normal groups. Kaplan-Meier curves allowed us to determine the impact of hub genes on melanoma prognosis.
The flow chart of this study is shown in Figure 1. In the training cohort, 17 immune cells associated with survival were ascertained using univariate COX regression analysis and subsequently screened by LASSO (Figure 2A, B) to obtain the ICRS constructed from Activated CD8 T cells and Immature B cells. The ICRS was calculated on account of the formula: the abundance of Activated CD8 T cell × (-0.429966) + the abundance of Immature B cell × (-1.514149). The distribution of patients, state of survival, and expression abundance of Activated CD8 T cells and Immature B cells between the different ICRS are shown in Figure 2C. Kaplan-Meier curves indicated that patients in the high-ICRS group had great shorter survival times than those in the low-ICRS group (P < 0.05) (Figure 2D). Next, immune characteristics between the different ICRS groups were analyzed using ESTIMATE and CIBERSORT methods. We noted that the immune score, ESTIMATE score and stromal score were great higher in the low -ICRS group (Figure 2E–G), while the tumor purity was considerably higher in the high-ICRS group (Figure 2H). Furthermore, it is worth noting that in the results of immune cell infiltration, both T cell CD8 and B cells native were discovered to be more highly expressed in the low immune group (Figure 2I, J), which is consistent with our previous analysis (P < 0.001).
To further determine the accuracy of the ICRS, GSE65904 was employed as the verification cohort for analysis in this study. The abundance of 23 immune cells was evaluated by ssGSEA and modeled using the same formula as the training cohort. Next, melanoma patients were divided into different ICRS groups on the strength of the same cut-off of the ICRS as the training cohort. The distribution of patients, state of survival and the abundance of Activated CD8 T cell and Immature B cell expression in the different ICRS groups are shown (Figure 3A). Survival curves show a better prognosis for patients in the high ICRS group (P < 0.05) (Figure 3B). Next, we noted that immune scores, ESTIMATE scores, and stromal scores were significantly higher in the low-ICRS group of the verification cohort (Figure 3C–E). Meanwhile, the purity of tumors in the high subgroup significantly exceeded that of the ground subgroup. (Figure 3F). Furthermore, most immune cells were significantly more elevated in the low-ICRS group (Figure 3G, H). These results are consistent with the training cohort, showing that the ICRS model identifies immune subtypes of melanoma well and effectively predicts the prognosis of melanoma patients.
Analyses were executed in the training cohort to elucidate the functional differences between the different ICRS groups. First, GSVA was utilized to illuminate the differences in pathways between the two groups. The different ICRS groups are mainly enriched in immune-related pathways such as allograft rejection and antigen processing and presentation (Figure 4A). Next, differential analysis between different ICRS groups was identified using the Wilcoxon rank sum test, and 972 differential genes were obtained (Figure 4B). GO or KEGG enrichment analysis of 972 differential genes to identify the functions of these genes or the pathways involved in regulation. The results of GO enrichment analysis suggested that in terms of biological processes (BP), differential genes were mainly involved in processes such as leukocyte-mediated immunity and lymphocyte-mediated immunity, and in terms of cellular components (CC), differential genes were associated with T cell receptor complexes, MHC protein complex, and other components. In terms of molecular functions (MF), differential genes are mainly immunoglobulin receptor binding, cytokine activity, chemokine receptor binding and other procedures (Figure 4D). KEGG pathways are primarily enriched in Cell adhesion molecules and hematopoietic cell lineage (Figure 4C). These differential terms may be potentially responsible for prolonged survival time in low-ICRS groups.
We further screened these DEGs using Lasso and random forest classifiers to filter the hub genes with predictive values for melanoma. LASSO regression analysis screened 45 significant candidate hub genes (Figure 5A). In the random forest, the number of random seeds is set to 123,456. We chose 500 trees as the random forest model's parameters, representing the model's stable error (Figure 5B). Finally, we screened the 43 hub genes with the highest Gini index scores. Finally, we filtered the 43 genes with the highest Gini index scores, and we obtained 45 hub genes in the Lasso analysis. The intersection of five potential candidate genes in RF and LASSO was shown by a Venn diagram (Figure 5C), and five genes (CXCL13, FCRL3, NCF1C, PLA2G2D, and SLAMF6) were identified as hub signature genes for subsequent analysis. Next, correlations among the five hub genes were analyzed using correlation analysis (Figure 5D).
To better investigate the role of 5 hub genes in melanoma prognosis, we analyzed scRNA-sequencing data, which included 4645 cells. A total of 4612 cells and 22,289 genes eventually passed the quality control filter and were utilized for further data analysis. Unsupervised clustering analysis was performed using the whole transcriptome of all single cells after batch effect correction, and the results were visualized using TSNE for downscaling. By using the "singR" package to annotate the cells, we obtained a total of 8 clusters of cells: T cells, B cell, Neurons, Monocyte, Tissue stem cells, NK cell, Fibroblasts, Endothelial cells, and visualized the top 5 significant differential genes in each cluster by heatmap (Figure 6A, B). Finally, we analyzed the distribution of the five hub genes in each cell cluster. We showed that SLAMF6 and FCRL3 were mainly distributed in T cells, B cells, and NK cell subpopulations, NCF1C was told primarily on B cells and Monocyte, CXCL9 was distributed in T cells, while PLA2G2D was not expressed in almost any cells (Figure 6C–G).
To better comprehend the intercommunication of different cell types in melanoma, we analyzed the interactions between eight cell types using cellular communication analysis. And further extracted individual cells to observe the communication between this cell and other cells (Figure 7A), and it was found that cellular communication between T cells and B cells was more assertive (Figure 7B). Next, we analyzed the cellular communication involved in the core gene. The results found that the core gene CXCL13 is present in the CXCL pathway between T cells and b cells, where CXCL13 is the ligand and CXCR5 is the ligand, so we further investigated the communication role and cellular action type analysis between immune cells in the CXCL pathway (Figure 7C). It was found that B cells, Endothelial cells, monocytes and Fibroblasts can act as both ligand and receptor cells in this pathway, where Endothelial cells, monocytes, and Fibroblasts can send signals to other cells. In contrast, T cells, Neurons, Tissue stem cells, and NK cells, can only act as ligand cells (Figure 7D, E). These results elucidate the possible interactions between these cell types and the cellular pathways involved in core genes. These help us further investigate the combined role of cell types and hub genes in melanoma prognosis.
To investigate the expression characteristics of hub genes between normal samples and melanoma, we downloaded 812 normal skin samples from the GTEx database for joint analysis with 463 melanoma samples from the TCGA database to ascertain the expression differences of hub genes between the tumor and normal groups. The five hub genes, CXCL13, FCRL3, NCF1C, PLA2G2D, and NCF1C, were highly expressed in normal samples with significant significance (Figure 8K–O). Next, further survival curve analysis showed that the five hub genes, CXCL13, FCRL3, PLA2G2D, NCF1C, and SLAMF6, had longer survival times for high expression. Although the gene PLA2G2D was insignificant in the verification group (Figure 8A–J). The above results suggest that these five genes are good prognostic genes and may play a key role in influencing the recurrence of melanoma.
Melanoma is composed of many different cells and the interaction between these cells is associated with tumor formation, metastasis, treatment resistance and immune infiltration [19]. These different cells interact with each other through ligand-receptor interactions. Because of the importance of this intercellular communication for patient prognosis, therapeutic approaches targeting intercellular interactions are becoming increasingly popular in clinical research. However, there have not been sufficient studies to identify the genes involved in intercellular communication.
In this study, the TCGA dataset as a training cohort, the abundance of 23 immune cells was assessed by ssGSEA, and a combination of univariate and LASSO COX regression analysis was utilized to screen for Activated CD8 T cells and Immature B cells, which were significantly relevant with melanoma prognosis. A validated ICRS was constructed on account of these two immune cells. We investigated the model's validity in predicting melanoma survival, showing that patients in the low-ICRS group had a better prognosis. The ICRS model also had an excellent diagnostic performance in the external verification cohort. These results indicate that prognostic models constructed on account of immune cells can precisely predict the survival rate of melanoma patients.
In both cell types, the risk factors for T and B cells are less than 0, indicating that these two immune cells are protective factors in melanoma. Alterations in immune cell content play a significant impact in melanoma progression. Studies have reported that B cell in patients with advanced melanoma is pivotal for an efficient anti-tumor immune response and stable disease control. At the same time, B cells could play a crucial role in regulating early tumor progression in melanoma patients [20]. On the other hand, activated CD8 T cells were found in our study to play a crucial role in melanoma development. An increased number of CD8+ T cells can control tumor growth [21].
With the development of bioinformatics, many risk models have been constructed to assess the prognostic outcome of melanoma patients. Song et al. created a predictive model that utilizes coagulation-related genes to predict patient survival and facilitate clinical treatment management [22]. Yang et al. reported that the analysis of cuproptosis-related lncRNAs in SKCM opens new directions for SKCM treatment [23]. Fu et al. constructed an immunogenic cell death(ICD) gene-related model to initially explore the ICD and possible signaling pathways related to melanoma prognosis [24]. Nonetheless, these studies have all been multi-gene combination analyses to predict the prognostic impact on melanoma. Still, constructing predictive models with two immune cells, T cells and B cells, to expect melanoma survival is not reported. Notably, a predictive model built with T and B cells could represent the immune profile of melanoma and afford a basis for assessing the response to clinical treatment of melanoma.
In addition, the difference genes between the different ICRS groups were obtained by The Wilcoxon rank sum test. The different genes were discovered to be mainly enriched in immune-related pathways by GO and KEGG enrichment analysis. GSVA also elucidated the pathway differences between the different ICRS groups. On account of this, we further screened the hub genes by Lasso and random forest algorithms, and five hub genes were screened: CXCL13, FCRL3, NCF1C, PLA2G2D and SLAMF6. In a previous study, CXCL13 was found to be a prognostic biomarker in melanoma patients treated with immune checkpoint inhibition [25], which enrolls B cells into the tumor microenvironment, and CXCL13 decreases regulatory B cells in the tumor microenvironment and restrains tumor growth [26]. CXCL13 was identified as an immune-related biomarker relevant to tumor formation and prognosis in patients with melanoma [27]. These are consistent with our findings. FCRL3 accelerates TLR9-induced B-cell activation and inhibits plasma cell decomposition [28] and has also been validated as an underlying prognostic biomarker relevant to melanoma metastasis [29]. NCF1 has been reported to cause an autosomal recessive disorder, chronic sarcoidosis [30]. Its association with cancer and melanoma has not been reported and deserves further exploration in the future. PLA2G2D is associated with melanoma metastasis, and its more in-depth mechanisms need to be further explored. Anti-SLAMF6 was found to be effective in reclaiming CD8+ T cell function and directly influencing tumor progression [31]. SLAMF6 and SLAM-related protein levels in CD8+ T lymphocytes were lower in patients with severe repletion anemia than in normal controls and were inversely associated with CD8+ T lymphocyte functional molecular levels [32]. SLAMF6 is a suppressive immune receptor, and its deletion allows powerful CD8+ T cells to eliminate tumors [33]. Its relationship with T cells is closely related, but the mechanism of its effect on melanoma has not been explored.
To further investigate the hub genes relevant to melanoma prognosis, we integrated bulk and scRNA sequencing data for analysis. In this report, the data set GSE72056 was obtained from the GEO database, including a total of 4645 cells from 19 melanoma patients [34]. The cells were annotated by singR software package, and we obtained eight cell populations of T cells, B cells, Neurons, Monocyte, tissue stem cells and NK cells. The present study used the interaction between T cells and B cells to construct the model. Tumor-associated B cells can attract effector T cells into the tumor and protect them. They can help activate the immune role of T cells, acting together with other immune components [35]. In addition, B cells may also activate and enroll other immune effector cells by secreting a range of cytokines (including TNF, IL-2, IL-6, and INF-γ) [36]. Several studies in cancer have now shown that the existence of TLS and B cells in tumor tissue is relevant to a better prognosis and that the anti-tumor effect of T cells is more substantial when B cells are present [37]. The tumor microenvironment also contains many immunosuppressive cells, such as mesenchymal cells, fibroblasts, regulatory T cells (Treg), macrophages (TAM), and myeloid-derived suppressor cells (MDSC). Treg cells can reduce anti-tumor immunity and accelerate tumor proliferation and metastasis by suppressing the functions of CD8+ T cells, NK cells, B cells and APC cells [38,39]. It has also been reported that melanoma growth and angiogenesis are precipitated after the inhibition of sensory neuron activity and decreased after the overstimulation of these neurons [40]. As a highly immunogenic tumor, the interrelationship between immune cells in melanoma has an essential impact on antitumor immune effects [41]. Therefore, studying the cellular composition of the tumor microenvironment and its dynamic changes is crucial for melanoma progression.
Then, we predicted potential ligand-receptor pairs between immune cells and found that the core genes screened by machine learning are mainly involved in the CXCL pathway during cellular communication. Analysis of the pathway's receptor-ligand pairs revealed that CXCL13-CXCR5 receptor-ligand pairs are only present in cellular communication between T and B cells. It was reported that CXCR5 is mainly conveyed in mature B cells and follicular helper T cells and that the CXCL13/CXCR5 axis coordinates cell-cell interactions and regulates lymphocyte infiltration in the tumor microenvironment [42]. Altered chemokine receptor-ligand expression of CXCR5/CXCL13 has been reported to be relevant in establishing B-cell dysfunction during HIV-1 infection [43]. The present study hypothesized that altered chemokine receptor-ligand expression of CXCR5/CXCL13 influences the immune infiltration of T lymphocytes during melanoma formation, which is one of the highlights of this study and requires further verification by subsequent studies. The present study also has some limitations. First, although we provided a prognosis model to predict melanoma survival, more prospective studies are needed to confirm the model's reliability. We also need more experimental evidence to prove our conclusions and elucidate the exact mechanisms of hub genes in melanoma development.
In this study, a novel ICRS was built to effectively predict the immune subtypes of SKCM patients. A machine learning algorithm further screened five hub genes that may be involved in regulating the immune microenvironment. In addition, on account of the analysis of scRNA-seq data from melanoma samples, this study described the cell clusters in melanoma, leading to a deeper understanding of underlying intercellular interactions in the immune microenvironment of melanoma, and identified ligand-receptor genes that determine intercellular interactions, suggesting that they may be underlying therapeutic targets affecting the prognosis of melanoma.
This study was supported by the Guangdong Provincial Science and Technology Innovation Strategy Special Fund (Climbing Plan Special Fund) Project (grant No. pdjh2020a0300), Key Laboratory of Pharmacovigilance Technology Research and Evaluation of State Drug Administration and Innovation Team of Health Economics and Social Security Research.
The authors declare that there is no conflict of interest.
| [1] |
T. Yamauchi, S. Shangraw, Z. Zhai, D. Ravindran Menon, N. Batta, R. P. Dellavalle, et al., Alcohol as a non-UV social-environmental risk factor for melanoma, Cancers (Basel), 14 (2022), 5010. https://doi:10.3390/cancers14205010 doi: 10.3390/cancers14205010
|
| [2] |
F. S. Hodi, S. J. O'Day, D. F. McDermott, R. W. Weber, J. A. Sosman, J. B. Haanen, et al., Improved survival with ipilimumab in patients with metastatic melanoma, N. Engl. J. Med., 363 (2010), 711–723. https://doi:10.1056/NEJMoa1003466 doi: 10.1056/NEJMoa1003466
|
| [3] |
Y. Qu, S. Zhang, Y. Zhang, X. Feng, F. Wang, Identification of immune-related genes with prognostic significance in the microenvironment of cutaneous melanoma, Virchows Arch., 478 (2021), 943–959. https://doi:10.1007/s00428-020-02948-9 doi: 10.1007/s00428-020-02948-9
|
| [4] |
Y. Xie, J. Zhang, M. Li, Y. Zhang, Q. Li, Y. Zheng, et al., Identification of lactate-related gene signature for prediction of progression and immunotherapeutic response in skin cutaneous melanoma, Front. Oncol., 12 (2022), 818868. https://doi:10.3389/fonc.2022.818868 doi: 10.3389/fonc.2022.818868
|
| [5] |
N. A. N. Jorge, J. G. V. Cruz, M. A. M. Pretti, M. H. Bonamino, P. A. Possik, M. Boroni, Poor clinical outcome in metastatic melanoma is associated with a microRNA-modulated immunosuppressive tumor microenvironment, J. Transl. Med., 18 (2020), 56. https://doi:10.1186/s12967-020-02235-w doi: 10.1186/s12967-020-02235-w
|
| [6] |
L. Fattore, C. F. Ruggiero, D. Liguoro, R. Mancini, G. Ciliberto, Single cell analysis to dissect molecular heterogeneity and disease evolution in metastatic melanoma, Cell Death Dis., 10 (2019), 827. https://doi:10.1038/s41419-019-2048-5 doi: 10.1038/s41419-019-2048-5
|
| [7] |
A. J. Combes, B. Samad, J. Tsui, N. W. Chew, P. Yan, G. C. Reeder, et al., Discovering dominant tumor immune archetypes in a pan-cancer census, Cell, 185 (2022), 184–203. https://doi:10.1016/j.cell.2021.12.004 doi: 10.1016/j.cell.2021.12.004
|
| [8] |
S. S. Potter, Single-cell RNA sequencing for the study of development, physiology and disease, Nat. Rev. Nephrol., 14 (2018), 479–492. https://doi:10.1038/s41581-018-0021-7 doi: 10.1038/s41581-018-0021-7
|
| [9] |
K. Asada, K. Takasawa, H. Machino, S. Takahashi, N. Shinkai, A. Bolatkan, et al., Single-cell analysis using machine learning techniques and its application to medical research, Biomedicines, 9 (2021), 1513. https://doi:10.3390/biomedicines9111513 doi: 10.3390/biomedicines9111513
|
| [10] |
M. P. Kumar, J. Du, G. Lagoudas, Y. Jiao, A. Sawyer, D. C. Drummond, et al., Analysis of single-cell RNA-Seq identifies cell-cell communication associated with tumor characteristics, Cell Rep., 25 (2018), 1458–1468. https://doi:10.1016/j.celrep.2018.10.047 doi: 10.1016/j.celrep.2018.10.047
|
| [11] |
J. Chen, S. Suo, P. P. Tam, J. J. Han, G. Peng, N. Jing, Spatial transcriptomic analysis of cryosectioned tissue samples with Geo-seq, Nat. Protoc., 12 (2017), 566–580. https://doi:10.1038/nprot.2017.003 doi: 10.1038/nprot.2017.003
|
| [12] |
L. Huang, C. Wu, D. Xu, Y. Cui, J. Tang, Screening of important factors in the early sepsis stage based on the evaluation of ssGSEA algorithm and ceRNA regulatory network, Evol. Bioinf. Online, 17 (2021), 11769343211058463. https://doi:10.1177/11769343211058463 doi: 10.1177/11769343211058463
|
| [13] |
Y. Li, J. Liu, X. Gao, B. Jie, M. Kim, P. T. Yap, et al., Multimodal hyper-connectivity of functional networks using functionally-weighted LASSO for MCI classification, Med. Image Anal., 52 (2019), 80–96. https://doi:10.1016/j.media.2018.11.006 doi: 10.1016/j.media.2018.11.006
|
| [14] |
K. Yoshihara, M. Shahmoradgoli, E. Martínez, R. Vegesna, H. Kim, W. Torres-Garcia, et al., Inferring tumour purity and stromal and immune cell admixture from expression data, Nat. Commun., 4 (2013), 2612. https://doi:10.1038/ncomms3612 doi: 10.1038/ncomms3612
|
| [15] |
X. Deng, T. Li, L. Mo, F. Wang, J. Ji, X. He, et al., Machine learning model for the prediction of prostate cancer in patients with low prostate-specific antigen levels: A multicenter retrospective analysis, Front. Oncol., 12 (2022), 985940. https://doi:10.3389/fonc.2022.985940 doi: 10.3389/fonc.2022.985940
|
| [16] |
Y. Zhou, W. Shi, D. Zhao, S. Xiao, K. Wang, J. Wang, Identification of immune-associated genes in diagnosing aortic valve calcification with metabolic syndrome by integrated bioinformatics analysis and machine learning, Front. Immunol., 13 (2022), 937886. https://doi:10.3389/fimmu.2022.937886 doi: 10.3389/fimmu.2022.937886
|
| [17] |
J. Vallejo, R. Saigusa, R. Gulati, S. S. Armstrong Suthahar, V. Suryawanshi, A. Alimadadi, et al., Combined protein and transcript single-cell RNA sequencing in human peripheral blood mononuclear cells, BMC Biol., 20 (2022), 193. https://doi:10.1186/s12915-022-01382-4 doi: 10.1186/s12915-022-01382-4
|
| [18] |
M. Efremova, M. Vento-Tormo, S. A. Teichmann, R. Vento-Tormo, CellPhoneDB: inferring cell-cell communication from combined expression of multi-subunit ligand-receptor complexes, Nat. Protoc., 15 (2020), 1484–1506. https://doi:10.1038/s41596-020-0292-x doi: 10.1038/s41596-020-0292-x
|
| [19] |
D. Hanahan, R. A. Weinberg, Hallmarks of cancer: the next generation, Cell, 144 (2011), 646–674. https://doi:10.1016/j.cell.2011.02.013 doi: 10.1016/j.cell.2011.02.013
|
| [20] |
A. D. Therien, G. M. Beasley, K. E. Rhodin, N. E. Farrow, D. S. Tyler, D. Boczkowski, et al., Spatial biology analysis reveals B cell follicles in secondary lymphoid structures may regulate anti-tumor responses at initial melanoma diagnosis, Front. Immunol., 13 (2022), 952220. https://doi:10.3389/fimmu.2022.952220 doi: 10.3389/fimmu.2022.952220
|
| [21] |
M. S. Ryu, M. Y. Woo, D. Kwon, A. E. Hong, K. Y. Song, S. Park, et al., Accumulation of cytolytic CD8+ T cells in B16-melanoma and proliferation of mature T cells in TIS21-knockout mice after T cell receptor stimulation, Exp. Cell Res., 327 (2014), 209–221. https://doi:10.1016/j.yexcr.2014.07.028 doi: 10.1016/j.yexcr.2014.07.028
|
| [22] |
B. Song, H. Chi, G. Peng, Y. Song, Z. Cui, Y. Zhu, et al., Characterization of coagulation-related gene signature to predict prognosis and tumor immune microenvironment in skin cutaneous melanoma, Front. Oncol., 12 (2022), 975255. https://doi:10.3389/fonc.2022.975255 doi: 10.3389/fonc.2022.975255
|
| [23] |
X. Yang, X. Wang, X. Sun, M. Xiao, L. Fan, Y. Su, et al., Construction of five cuproptosis-related lncRNA signature for predicting prognosis and immune activity in skin cutaneous melanoma, Front. Genet., 13 (2022), 972899. https://doi:10.3389/fgene.2022.972899 doi: 10.3389/fgene.2022.972899
|
| [24] |
W. Fu, G. Ma, Significance of immunogenic cell death-related genes in prognosis prediction and immune microenvironment landscape of patients with cutaneous melanoma, Front. Genet., 13 (2022), 988821. https://doi:10.3389/fgene.2022.988821 doi: 10.3389/fgene.2022.988821
|
| [25] |
C. S. Groeneveld, J. Fontugne, L. Cabel, I. Bernard-Pierrot, F. Radvanyi, Y. Allory, et al., Tertiary lymphoid structures marker CXCL13 is associated with better survival for patients with advanced-stage bladder cancer treated with immunotherapy, Eur. J. Cancer, 148 (2021), 181–189. https://doi:10.1016/j.ejca.2021.01.036 doi: 10.1016/j.ejca.2021.01.036
|
| [26] |
L. Shen, J. Li, Q. Liu, M. Das, W. Song, X. Zhang, et al., Nano-trapping CXCL13 reduces regulatory B cells in tumor microenvironment and inhibits tumor growth, J. Control. Release, 343 (2022), 303–313. https://doi:10.1016/j.jconrel.2022.01.039 doi: 10.1016/j.jconrel.2022.01.039
|
| [27] |
Z. Si, H. Hu, Identification of CXCL13 as an Immune-related biomarker associated with tumorigenesis and prognosis in cutaneous melanoma patients, Med. Sci. Monit., 27 (2021), e932052. https://doi:10.12659/msm.932052 doi: 10.12659/msm.932052
|
| [28] |
F. J. Li, D. M. Schreeder, R. Li, J. Wu, R. S. Davis, FCRL3 promotes TLR9-induced B-cell activation and suppresses plasma cell differentiation, Eur. J. Immunol., 43 (2013), 2980–2992. https://doi:10.1002/eji.201243068 doi: 10.1002/eji.201243068
|
| [29] |
Y. Li, S. Lyu, Z. Gao, W. Zha, P. Wang, Y. Shan, et al., Identification of potential prognostic biomarkers associated with cancerometastasis in skin cutaneous melanoma, Front. Genet., 12 (2021), 687979. https://doi:10.3389/fgene.2021.687979 doi: 10.3389/fgene.2021.687979
|
| [30] |
D. B. Kuhns, A. P. Hsu, D. Sun, K. Lau, D. Fink, P. Griffith, et al., NCF1 (p47(phox))-deficient chronic granulomatous disease: comprehensive genetic and flow cytometric analysis, Blood Adv., 3 (2019), 136–147. https://doi:10.1182/bloodadvances.2018023184 doi: 10.1182/bloodadvances.2018023184
|
| [31] |
B. Yigit, N. Wang, E. T. Hacken, S. S. Chen, A. K. Bhan, A. Suarez-Fueyo, et al., SLAMF6 as a regulator of exhausted CD8(+) T cells in cancer, Cancer Immunol. Res., 7 (2019), 1485–1496. https://doi:10.1158/2326-6066.Cir-18-0664 doi: 10.1158/2326-6066.Cir-18-0664
|
| [32] |
B. Liu, L. Zeng, Y. Shao, R. Fu, Expression and function of SLAMF6 in CD8(+) T lymphocytes of patients with severe aplastic anemia, Cell Immunol., 364 (2021), 104343. https://doi:10.1016/j.cellimm.2021.104343 doi: 10.1016/j.cellimm.2021.104343
|
| [33] |
E. Hajaj, G. Eisenberg, S. Klein, S. Frankenburg, S. Merims, I. B. David, et al., SLAMF6 deficiency augments tumor killing and skews toward an effector phenotype revealing it as a novel T cell checkpoint, Elife, 9 (2020). https://doi:10.7554/eLife.52539 doi: 10.7554/eLife.52539
|
| [34] |
I. Tirosh, B. Izar, S. M. Prakadan, M. H. Wadsworth, D. Treacy, J. J. Trombetta, et al., Dissecting the multicellular ecosystem of metastatic melanoma by single-cell RNA-seq, Science, 352 (2016), 189–196. https://doi:10.1126/science.aad0501 doi: 10.1126/science.aad0501
|
| [35] |
J. Griss, W. Bauer, C. Wagner, M. Simon, M. Chen, K. Grabmeier-Pfistershammer, et al., B cells sustain inflammation and predict response to immune checkpoint blockade in human melanoma, Nat. Commun., 10 (2019), 4186. https://doi:10.1038/s41467-019-12160-2 doi: 10.1038/s41467-019-12160-2
|
| [36] |
W. H. Fridman, F. Petitprez, M. Meylan, T. W. Chen, C. M. Sun, L. T. Roumenina, et al., B cells and cancer: To B or not to B, J. Exp. Med., 218 (2021), e20200851. https://doi:10.1084/jem.20200851 doi: 10.1084/jem.20200851
|
| [37] |
M. Lauss, M. Donia, I. M. Svane, G. Jönsson, B cells and tertiary lymphoid structures: Friends or foes in cancer immunotherapy, Clin. Cancer Res., 28 (2022), 1751–1758. https://doi:10.1158/1078-0432.Ccr-21-1130 doi: 10.1158/1078-0432.Ccr-21-1130
|
| [38] |
T. Gambichler, M. Bindsteiner, S. Höxtermann, S. Terras, A. Kreuter, Circulating CD4+ CD25 (high) CD127 (low) regulatory T cells are an independent predictor of advanced melanoma, Pigment Cell Melanoma Res., 26 (2013), 280–283. https://doi:10.1111/pcmr.12055 doi: 10.1111/pcmr.12055
|
| [39] |
C. Miracco, V. Mourmouras, M. Biagioli, P. Rubegni, S. Mannucci, I. Monciatti, et al., Utility of tumour-infiltrating CD25+FOXP3+ regulatory T cell evaluation in predicting local recurrence in vertical growth phase cutaneous melanoma, Oncol. Rep., 18 (2007), 1115–1122. https://doi.org/10.3892/or.18.5.1115 doi: 10.3892/or.18.5.1115
|
| [40] |
P. A. C. Costa, W. N. Silva, P. Prazeres, C. C. Picoli, G. D. A. Guardia, A. C. Costa, et al., Chemogenetic modulation of sensory neurons reveals their regulating role in melanoma progression, Acta Neuropathol. Commun., 9 (2021), 183. https://doi:10.1186/s40478-021-01273-9 doi: 10.1186/s40478-021-01273-9
|
| [41] |
S. Kalaora, A. Nagler, J. A. Wargo, Y. Samuels, Mechanisms of immune activation and regulation: lessons from melanoma, Nat. Rev. Cancer, 22 (2022), 195–207. https://doi:10.1038/s41568-022-00442-9 doi: 10.1038/s41568-022-00442-9
|
| [42] |
M. Hussain, D. Adah, M. Tariq, Y. Lu, J. Zhang, J. Liu, CXCL13/CXCR5 signaling axis in cancer, Life Sci., 227 (2019), 175–186. https://doi:10.1016/j.lfs.2019.04.053 doi: 10.1016/j.lfs.2019.04.053
|
| [43] |
A. Cagigi, F. Mowafi, L. V. Phuong Dang, K. Tenner-Racz, A. Atlas, S. Grutzmeier, et al., Altered expression of the receptor-ligand pair CXCR5/CXCL13 in B cells during chronic HIV-1 infection, Blood, 112 (2008), 4401–4410. https://doi.org/10.1182/blood-2008-02-140426 doi: 10.1182/blood-2008-02-140426
|
Figures(8)
Linqian Guo, Qingrong Meng, Wenqi Lin, Kaiyuan Weng. Identification of immune subtypes of melanoma based on single-cell and bulk RNA sequencing data[J]. Mathematical Biosciences and Engineering, 2023, 20(2): 2920-2936. doi: 10.3934/mbe.2023138
DownLoad:
