Single-cell RNA-seq profiling, clustering, annotation and identifications of macrophage markers
The flow chart of this study could be seen in Fig. 1. After processing scRNA-seq data based on the strict quality control metrics as previously mentioned, we obtained 16,291 cell samples with superior quality isolated from 48 melanoma tissues in GSE120575. A strong positive correlation between detected gene numbers and sequencing depth was observed with the Pearson's correlation of 0.97 (Fig. 2A). After normalizing the data, we selected the top 2,000 highly variable genes for the subsequent analysis (Fig. 2B). We conducted PCA using the top 2,000 variable genes to reduce the dimensionality, and the top 15 principal components with p value < 0.05 were selected for follow-up analysis. We subsequently adopted UMAP technique on the 15 subpopulations to visualize the high dimensional scRNA-seq data, and successfully identified 12 subclasses (Fig. 2C), which were later annotated to five cell types, including CD8 + T cells (CD8A + CD8B + CD3D + CD3E+), NK/NKT cells (GNLY + NKG7 + KLRD1+), CD4 + T cells (CD4 + CD3D + CD3E+), B cells (CD19 + MS4A1) and macrophages (CSF1R + CD68+), according to well-recognized marker genes (Fig. 2D, Figure S1). Macrophages were also found to have distinct gene expression profiles, with 865 genes differentially expressed according to |log2 FC|> 0.585 and adjusted p value < 0.05, which were identified as melanoma-related macrophage marker genes. The most common markers of each cell type are visualized as a Fig. 2E. In addition, the macrophage cluster was manually verified through well-acknowledged macrophage markers including FCGR3A, ITGAM, LST1, CSF1R, CD14, LYZ, TLR2, TLR4, S100A8, S100A9, CD80, CD86, CD68, CD163 (17–24) (Fig. 3).
GO, KEGG analyses and consensus analysis of macrophage marker genes
We conducted GO and KEGG enrichment analyses to investigate the correlative functions and pathways for the above 865 macrophage marker genes. As shown in Fig. 4A, these genes were mostly enriched in immunity-related features of macrophages, such as lysosomal membrane, lytic vacuole membrane, phagocytosis, positive regulation of cytokine production, regulation of cell-cell adhesion, leukocyte cell-cell adhesion, mononuclear cell differentiation, regulation of T cell activation, leukocyte mediated immunity, positive regulation of leukocyte activation, regulation of leukocyte cell-cell adhesion, activation of immune response, positive regulation of lymphocyte activation, lymphocyte differentiation. On the other hand, immunity-related pathways of macrophages, like phagosome, lysosome, Th1 and Th2 cell differentiation, cytokine-cytokine receptor interaction, Th17 cell differentiation, NF-kappa B, T cell receptor, PD-L1 expression and PD-1 checkpoint, chemokine, hematopoietic cell lineage, human immunodeficiency virus 1 infection, natural killer cell mediated cytotoxicity, cell adhesion molecules, neutrophil extracellular trap formation, etc. were significantly correlated with these genes (Fig. 4B).
Based on the above melanoma-related macrophage marker genes, we performed consistent clustering to cluster TCGA-SKCM samples and obtained 6 macrophage subtypes (Fig. 4C). Among the 6 macrophage subtypes, there were the most significantly different in survival and macrophage marker gene expressed between macrophage subtypes A and B. Specifically, macrophage subtype B with the highest expression of macrophage marker genes was associated with the optimal prognosis. On the contrary, macrophage subtype A with the lowest expression of macrophage marker genes was related to the worst prognosis (Fig. 4D, 4E).
We found that the expression of antigen-presenting related genes (including HLA-A, HLA-C, HLA-DPA1, HLA-DPB1, HLA-DPB2, HLA-DQA1, HLA-DQA2, HLA-DQB1, HLA-DQB2, HLA-DRA, HLA-DRB1, HLA-DRB5, HLA-DRB6), cytokine (involving IL2, IL4, IL5, IL6, IL10, IL12A, IL12B, IL13, IL15, IL18), chemokine (like CCL1, CCL2, CCL5, CCL17, CCL18, CCL22, CCL24, CXCL9, CXCL10, CXCL11, CXCL13, CXCL16), interferon family genes (such as IFNB1, IFNA1) and tumor necrosis factor were significantly different among the six macrophage subtypes, which meant that these subtypes of macrophages especially subtypes A and B where the differences were most striking, had significant biological behavioral differences in the anti-tumor immune process (Fig. 5).
In addition, we further analyzed the pathway enrichment characteristics of these macrophage subpopulations, and we discovered that anti-tumor immune-related pathways like antigen processing and presentation, natural killer cell mediated cytotoxicity, toll like receptor signaling pathway, nod like receptor signaling pathway, T cell receptor signaling pathway, B cell receptor signaling pathway, cytokine-cytokine receptor interaction, chemokine signaling pathway, focal adhesion, and classical tumor-related pathways such as MAPK signaling pathway, JAK-STAT signaling pathway were most remarkably enriched in B and C subtypes, followed by D and F subtypes, E subtype, and A subtype. Moreover, F subtype had the most pathway enrichment in lysosome and endocytosis, while the least enrichment in ubiquitin mediated proteolysis, regulation of autophagy, apoptosis. Thus, in-depth exploring the molecular differences between subtypes A and B was particularly necessary since immune- and tumor-related pathways were the most significantly enriched in A subtype, while the least involved in B subtype, suggesting that subgroup A was active, and yet subgroup B was silent in anti-tumor immune process (Fig. 6).
Establishment the prognostic signature
Then, 476 statistically differentially expressed genes were identified from macrophage subtypes A and B. To construct a macrophage-based immune-related risk signature, we first treated the TCGA-SKCM cohort as the training set, conducted univariate Cox regression analysis and K-M analysis on the 476 genes, and preliminarily screened out 327 genes with survival differences. Meanwhile, we downloaded a list of immune genes, including 2483 immune genes, from the IMMPORT website, intersected the immune genes with the 327 prognostic genes, ultimately picked out 102 immune related prognostic genes, which were involved in constructing a risk signature with stepwise multivariate Cox regression analysis (Fig. 7A, 7B). The final 21-gene signature formula was as follows: Risk score = SLA*0.668 + TNFRSF17*0.125 + NOD2*0.201 + CCL8*(-0.152) + CD1B*(-0.243) + ICOS*0.284 + PTPRC*(-0.321) + PIK3CG*(-0.176) + CD247*0.380 + IL18RAP*(-0.173) + CXCR6*(-0.255) + PTGER4*0.200 + CD72*(-0.407) + PRF1*(-0.254) + CXCL11*0.144 + HLA-DQA1*(-0.223) + KLRC1*(-0.178) + LCK*0.356 + HLA-DOB*(-0.250) + IL2RA*(-0.186) + IL15*(-0.167) (Fig. 7C).
Examination of the prognostic signature
Then, we calculated the risk score for each sample in the TCGA cohort and divided them into the high- (n = 220) and low-risk (n = 234) groups based on the optimal cut-off value. Figure 7D exhibited the characteristics of risk scores, survival status and signature gene expression of training set. Kaplan-Meier analysis revealed that low-risk patients had significantly better survival outcomes compared with high-risk patients (Fig. 7E). The 5-year areas under the curve (AUC) value for the TCGA-SKCM cohort was 0.691 (Fig. 7F). Moreover, the patients in the low-risk group also had better OS in different clinical subgroups (age (≤ 65 years old), sex (male versus female), T stage (T1-2 versus T3-4), N stage (N0 versus N1-3), and pathological stage (stage I-II versus stage III-IV)) (Figure S2).
Apparently, similar distributions of survival status, risk scores, and gene expression feature were observed in external validation sets GSE65904 and GSE59455, which proved that our risk model was of great prognostic prediction power and widely applicable to melanoma patients. Survival curves in the validation cohort GSE65904 and GSE59455 were also consistent with training cohort. Furthermore, the AUCs for 5-year OS to be 0.712 for GSE65904 and 0.697 for GSE59455, respectively. These results supported that our constructed risk signature had moderate potential for monitoring survival (Fig. 8A-8F).
Meta-analysis of the average risk score for the three cohorts showed significant heterogeneity (I2 = 95%, P < 0.05). While all high-risk groups had higher average risk score than low-risk groups (MD: 1.42, 95%CI: 1.21–1.62 in TCGA-SKCM; MD: 2.06, 95%CI: 1.72–2.40 in GSE65904; MD: 0.42, 95%CI: 0.06–0.78 in GSE59455). Overall effected test of the average risk score for three cohorts indicated a higher average risk score in high-risk group compared to low-risk group (P = 0.001) (Fig. 8G). There was no significant heterogeneity about the five-year survival events between high- and low-risk groups among the three cohorts (I2 = 48%, P = 0.14), and the 5-year survival rate of the low-risk groups was significantly higher than that of the high-risk groups (RR: 1.59, 95%CI: 1.22–2.08 in TCGA-SKCM; RR: 2.30, 95%CI: 1.15–4.59 in GSE65904; RR: 3.36, 95%CI: 1.55–7.27 in GSE59455). Overall effected test of the five-year survival events revealed more 5-year survival events in low-risk group than high-risk group (P < 0.00001) (Fig. 8H). These results demonstrated that our risk model was a reliable prognostic indicator for patients with melanoma.
Validation of the signature genes
Not only did these 21 genes influence the prognosis of melanoma patients (Figure S3), but were also closely related to the expression of PD-L1 gene (R: 0.49–0.79, p < 0.05) (Figure S4) and macrophage marker genes like CD80, CD86, CD163, CD204, CD206, CD209, CD301 (R: 0.37–0.94, p < 0.05). Besides, these genes also interacted closely with each other (R: 0.38–0.96, p < 0.05) to affect the occurrence and development of melanoma, and anti-tumour immune process (Figure S5). Further exploring the distribution of these signature genes in our annotated immune cells and we found that SLA, PTPRC, PTGER4, HLA-DQA1, CD72, PIK3CG were mostly expressed on macrophages and many other immune cells. ICOS, CD247, CXCR6, LCK, and PRF1 were mainly expressed in T cells and NK/NKT cells (Figure S6).
For melanoma patients treated with anti-PD1, anti-CTLA4 or anti-PD1 combined with anti-CTLA4 therapy, SLA, ICOS, PTPRC, PIK3CG, CD247, IL18RAP, CXCR6, CD72, PRF1, HLA-DQA1, LCK and HLA-DOB had higher gene expression in response groups than non-response groups. When melanoma patients treated with anti-PD1 or anti-PD1 combined with anti-CTLA4 therapy, NOD2, PTGER4 and IL2RA also had higher expression in response groups than non-response groups. On the other hand, for melanoma patients treated with anti-PD1 therapy, TNFRSF17 and CCL8 had higher gene expression in non-response groups than response groups. KLRC1 had higher gene expression in non-response groups than response groups for melanoma patients treated with anti-PD1 combined with anti-CTLA4 therapy as well. And CXCL11 also had higher gene expression in non-response groups than response groups for melanoma patients treated with anti-PD1, anti-CTLA4 or anti-PD1 combined with anti-CTLA4 therapy (Table 1).
Table 1
Signature gene expression differences between immunotherapy response and non-response melanoma patients.
Gene
|
Study
|
Anti Target
|
Response Mean
|
Non-Response Mean
|
Log2FC
|
FDR
|
P- Value
|
SLA
|
ERP105482,SRP011540,SRP070710, SRP094781,SRP150548,SRP230414, SRP250849,SRP302761
|
anti-PD1/anti-CTLA4/anti-PD1 + anti-CTLA4
|
85.5
|
64
|
0.732
|
0.001
|
0
|
TNFRSF17
|
SRP150548
|
anti-PD1
|
0
|
12.5
|
-5.733
|
0.324
|
0.024
|
NOD2
|
ERP105482
|
anti-PD1 + anti-CTLA4
|
773
|
335
|
1.234
|
0.077
|
0.001
|
GSE93157
|
anti-PD1
|
6.91
|
6.06
|
0.848
|
0.323
|
0.025
|
GSE122220
|
anti-PD1
|
6.64
|
5.21
|
1.43
|
0.992
|
0.049
|
CCL8
|
ERP105482
|
anti-PD1 + anti-CTLA4
|
74
|
26
|
2.293
|
0.014
|
0
|
SRP150548
|
anti-PD1
|
167
|
71.5
|
2.042
|
0.029
|
0
|
SRP070710
|
anti-PD1
|
105
|
302.5
|
-1.881
|
0.057
|
0.001
|
SRP230414
|
anti-PD1
|
99.5
|
321
|
-1.747
|
0.122
|
0.002
|
ICOS
|
ERP105482,SRP011540,SRP070710, SRP094781,SRP150548,SRP230414, SRP250849,SRP302761
|
anti-PD1/anti-CTLA4/anti-PD1 + anti-CTLA4
|
153
|
84.5
|
0.642
|
0.017
|
0.001
|
PTPRC
|
ERP105482,SRP011540,SRP070710, SRP094781,SRP150548,SRP230414, SRP250849,SRP302761
|
anti-PD1/anti-CTLA4/anti-PD1 + anti-CTLA4
|
9,618.50
|
5,679.00
|
0.553
|
0.015
|
0.001
|
PIK3CG
|
ERP105482,SRP011540,SRP070710, SRP094781,SRP150548,SRP230414, SRP250849,SRP302761
|
anti-PD1/anti-CTLA4/anti-PD1 + anti-CTLA4
|
718
|
451
|
0.436
|
0.042
|
0.004
|
CD247
|
ERP105482,SRP011540,SRP070710, SRP094781,SRP150548,SRP230414, SRP250849,SRP302761
|
anti-PD1/anti-CTLA4/anti-PD1 + anti-CTLA4
|
252
|
164
|
0.922
|
0
|
0
|
IL18RAP
|
ERP105482,SRP011540,SRP070710, SRP094781,SRP150548,SRP230414, SRP250849,SRP302761
|
anti-PD1/anti-CTLA4/anti-PD1 + anti-CTLA4
|
46.5
|
35
|
0.385
|
0.176
|
0.04
|
CXCR6
|
ERP105482,SRP011540,SRP070710, SRP094781,SRP150548,SRP230414, SRP250849,SRP302761
|
anti-PD1/anti-CTLA4/anti-PD1 + anti-CTLA4
|
210.5
|
187
|
0.415
|
0.055
|
0.006
|
PTGER4
|
SRP011540
|
anti-PD1
|
260
|
217.5
|
0.523
|
0.21
|
0.033
|
SRP250849
|
anti-PD1
|
65
|
35.33
|
1.426
|
0.737
|
0.037
|
CD72
|
ERP105482,SRP011540,SRP070710, SRP094781,SRP150548,SRP230414, SRP250849,SRP302761
|
anti-PD1/anti-CTLA4/anti-PD1 + anti-CTLA4
|
204
|
196
|
0.558
|
0.013
|
0.001
|
PRF1
|
ERP105482,SRP011540,SRP070710, SRP094781,SRP150548,SRP230414, SRP250849,SRP302761
|
anti-PD1/anti-CTLA4/anti-PD1 + anti-CTLA4
|
314.5
|
157
|
0.858
|
0.001
|
0
|
GSE93157
|
anti-PD1
|
8.27
|
7.08
|
1.19
|
0.323
|
0.034
|
CXCL11
|
ERP105482,SRP011540,SRP070710, SRP094781,SRP150548,SRP230414, SRP250849,SRP302761
|
anti-PD1/anti-CTLA4/anti-PD1 + anti-CTLA4
|
77.5
|
79
|
0.726
|
0.018
|
0.001
|
HLA-DQA1
|
ERP105482,SRP011540,SRP070710, SRP094781,SRP150548,SRP230414, SRP250849,SRP302761
|
anti-PD1/anti-CTLA4/anti-PD1 + anti-CTLA4
|
7,904.00
|
5,871.50
|
0.453
|
0.073
|
0.009
|
GSE93157
|
anti-PD1
|
9.53
|
5.9
|
3.634
|
0.323
|
0.02
|
KLRC1
|
ERP105482
|
anti-PD1
|
94
|
23.5
|
1.775
|
0.061
|
0.002
|
ERP105482
|
anti-PD1 + anti-CTLA4
|
78
|
115
|
1.253
|
0.424
|
0.049
|
LCK
|
ERP105482,SRP011540,SRP070710, SRP094781,SRP150548,SRP230414, SRP250849,SRP302761
|
anti-PD1/anti-CTLA4/anti-PD1 + anti-CTLA4
|
274.5
|
171
|
1.021
|
0
|
0
|
HLA-DOB
|
ERP105482,SRP011540,SRP070710, SRP094781,SRP150548,SRP230414,SRP250849,SRP302761
|
anti-PD1/anti-CTLA4/anti-PD1 + anti-CTLA4
|
30
|
25.5
|
0.879
|
0.003
|
0
|
IL2RA
|
SRP150548
|
anti-PD1
|
272
|
52.5
|
2.762
|
0.04
|
0
|
SRP094781
|
anti-PD1
|
456.5
|
276
|
1.198
|
0.314
|
0.014
|
Compared to the normal individuals, the expression of CCL8, CXCL11, CXCR6, PTPRC, PRF1, CD247, ICOS, LCK, HLA-DQA1, IL2RA were significantly overexpressed in melanoma patients. On the contrary, NOD2 and PTGER4 were highly expressed in normal individuals (Fig. 9). Additionally, we took advantage of immunohistochemistry results from the HPA database to further explore the protein expression levels of these signature genes in melanoma. Notably, CD247, HLA-DOB, IL18RAP, PIK3CG were observed mild upregulation and CXCL11, HLA-DQA1, ICOS, IL15, NOD2 moderate overexpression in melanoma tissues compared to that in non-tumour tissues. While PTPRC, SLA expression levels were mildly downregulated and CD1B was moderately downregulated in melanoma tissues. However, there was no significant difference in the expression levels of IL2RA, LCK, PTGER4, CD72 and PRF1 among which IL2RA, LCK were mildly expressed, PTGER4 was moderately expressed and CD72 and PRF1 were not expressed in both normal and melanoma tissues (Fig. 10). No data was available for TNFRSF17, CCL8, CXCR6, KLRC1 in the HPA database.
Tumour immune cell infiltration, tumour microenvironment, and immunotherapy response
Based on the macrophage related genes that were generated the 21-gene signature, low-risk group had pronouncedly higher expression of macrophage marker genes like TLR2, TLR4, CD80, CD86, CD14, CD163, CD209, CD68 and CSF1R in both TCGA-SKCM and GEE65904 cohorts (Fig. 11A). In terms for immune microenvironment, analogous results were noticed that low-risk groups had considerably higher levers of immune score, stromal score and estimate score in these two cohorts (Fig. 11B). Consistent findings were achieved that low-risk groups had obviously higher infiltration of M1 macrophages, activated memory CD4 + T cells, plasma cells, resting dendritic cells and γδ T cells, while memory B cells, M0 macrophages were higher infiltration in high-risk groups in these two cohorts. In addition, patients with low-risk scores had higher infiltration of CD8 + T cells and lower infiltration of resting memory CD4 + T cells and resting NK cells in TCGA-SKCM cohort. In GSE65904, patients with low-risk scores had lower infiltration of regulatory T cells (Fig. 11C). Moreover, CTLA4, PD-1 and PD-L1 expression were significantly higher in low-risk groups in these two cohorts (Fig. 11D). Correspondingly, low-risk group had a better response to anti-CTLA4 and anti-PD-1 therapy (Fig. 11E).
Functional and pathways assessment of risk signature
GSEA was used to explore the function and pathway enrichment of the high- and low-risk groups in the integrated cohorts of TCGA-SKCM and GSE65904. As shown in Fig. 12A, immune response to tumour cell, CD4 + αꞵ T cell differentiation, cytokine production involved in immune response, positive regulation of T cell migration, regulation of T cell mediated immunity, T cell mediated immunity, positive regulation of T cell proliferation, natural killer cell mediated immunity, coreceptor activity, abnormal B cell morphology, etc. biological activities were profoundly enriched in low-risk group. Likewise, low-risk group was strikingly enriched in immune-associated pathways, such as T cell receptor, primary immunodeficiency, antigen processing and presentation, natural killer cell mediated cytotoxicity, graft versus host disease, cell adhesion molecules cams, cytokine-cytokine receptor interaction, toll like receptor, nod like receptor signalling pathway etc. While rare immune functions and pathways were enriched in high-risk patients (Fig. 12).