Development and Validation of a CD8+ T Cell Infiltration-Related Signature for Melanoma Patients

Abstract

Aim: Immunotherapy shows efficacy in only a subset of melanoma patients. Here, we intended to construct a risk score model to predict melanoma patients' sensitivity to immunotherapy.

Methods: Integration analyses were performed on melanoma patients from high-dimensional public datasets. The CD8+ T cell infiltration related genes (TIRGs) were selected via TIMER and CIBERSORT algorithm. LASSO Cox regression was performed to screen for the crucial TIRGs. Single sample gene set enrichment analysis (ssGSEA) and ESTIMATE algorithm were used to evaluate the immune activity. The prognostic value of the risk score was determined by univariate and multivariate Cox regression analysis.

Results: 184 candidate TIRGs were identified in melanoma patients. Based on the candidate TIRGs, melanoma patients were classified into three clusters which were characterized by different immune activity. Six signature genes were further screened out of 184 TIRGs and a representative risk score for patient survival was constructed based on these six signature genes. The risk score served as an indicator for the level of CD8+ T cell infiltration and acted as an independent prognostic factor for the survival of melanoma patients. By using the risk score, we achieved a good predicting result for the response of cancer patients to immunotherapy. Moreover, pan-cancer analysis revealed the risk score could be used in a wide range of non-hematologic tumors.

Conclusions: Our results showed the potential of using signature gene-based risk score as an indicator to predict melanoma patients' sensitivity to immunotherapy.

Keywords: CD8+ T cells; immune response; immunotherapy; melanoma; single cell RNA sequencing analysis.

Figure 1

Identification and enrichment analysis of TIRGs in melanoma cells. (A) Venn diagram to identify TIRGs in melanoma patients from the TCGA dataset by using TIMER and CYBERSORT. (B) Venn diagram to screen TIRGs in melanoma patients from datasets TGCA, GSE65094 and GSE19234. (C) t-SNE analysis of single cell sequencing dataset GSE72056 illustrates gene expression patterns in different cell types (shown in different colors). (D, E) The expression pattern of CD8A (D) and CD8B (E) indicates the distribution of CD8+ T cells. (F–H) Illustration of three excluded TIRGs that were highly expressed in more than 50% of CD8+ T cells but less than 10% of the remaining cells in TME: CTSW (F), DTHD1 (G) and LAG3 (H). (I) KEGG enrichment analysis of the candidate TIRGs in melanoma TME.

Figure 2

Stratificaton of melanoma patients via the expression of 184 TIRGs. (A) 448 melanoma patients from TCGA dataset are classified into three clusters based on the selected TIRGs by optima selection of unsupervised clustering. (B) Principle component analysis on the expression level of 184 TIRGs. (C) Clinical characteristics and RNA expression level of 184 TIRGs of melanoma patients from cluster 1, 2 and 3. (D, E) ESTIMATE analysis of immune score (D) and stromal score (E) shows a significant difference among three clusters in melanoma patients. (F, G) Melanoma patients in cluster 2 and 3 have significantly longer overall survival (F) and longer disease-specific survival (G) than those in cluster 1.

Figure 3

Construction of a risk score based on six signature genes. (A, B) The LASSO Cox regression model was constructed from 184 signature genes, and the tuning parameter (λ) was calculated based on the partial likelihood deviance with ten-fold cross validation. The six signature genes were identified according to the best fit profile. (C) Melanoma patients in the TCGA training set were divided into two populations according to the median value of the risk scores. (D–F) Melanoma patients with lower risk score have significantly longer survival in the TCGA whole set (D), GSE65094 data sets (E) and GSE22153 data sets (F).

Figure 4

Characterization of risk score with immune activity. (A, B) Correlation between risk score and the infiltrating number of CD8 T cells in melanoma patients by analysis with TIMER (A) and CIBERSORT (B). (C) Comparison of seven T cell marker expression level between melanoma patients with high and low risk scores in the TCGA cohort. (D) Relationship between risk score and CD8+ T cell activation related genes. Adjusted P values were showed as: *P < 0.05; **P < 0.01; ***P < 0.001, ****P < 0.0001. (E) GSEA of gene expression from data sets TCGA_SKCM shows genes involving in CD8+ T cell infiltration signaling pathway mostly enriched in melanoma patients with low risk scores.

Figure 5

Analysis of mutation burden in high and low risk group. (A, B) Mutation landscape of 222 melanoma tumor samples with low risk (A) and high risk (B). Central matrix shows somatic mutations with colors indicating different types of mutations and genes mutated at high frequency are represented in the left list. The top bar plot shows the number of gene mutations in each sample and the mutation rate of significantly mutated genes is displayed on the right. (C) Comparison of mutation frequency between high risk group and low risk group. (D) Comparison of tumor mutation burden between high risk group and low risk group.

Figure 6

Integration analysis of risk score and clinical characteristics to predict the survival of melanoma patients. (A) Univariate analysis shows only gender is not significantly correlated with disease progression. (B) Multivariate analysis shows risk score, breslow, and age is significantly correlated with disease progression. (C) Nomogram including risk score constructed to predict the 1-, 3-, and 5-year survival of melanoma patients in the TCGA cohort. (D, E) Calibration curve of the nomogram for predicting the probability of OS at 1 and 3 years.

Figure 7

Risk score to predict the efficacy of immunotherapy on cancer patients. (A–E) Expression of immunotherapy targeted-genes in low risk and high risk melanoma patients from datasets TCGA_SKCM and GSE65904. (A) PDCD1-related genes. (B) CTLA4-related genes. (C) TIM3-related genes. (D) LAG3-related genes. (E) TIGIT-related genes. (F, G) Melanoma patients from data set GSE35640 who received immunotherapy but have no response (non-responder) show higher risk scores than those responders. (K) ROC curve showing the performance of our model for predicting the efficacy of immunotherapy on urothelial cancer patients in data set IMvigor210 at all classification thresholds (AUC=0.604) at all classification thresholds (AUC=0.753). (I, J) Non-responding urothelial cancer patients from data set IMvigor210 show higher risk scores than those responders. (H) ROC curve showing the performance of our model for predicting the efficacy of immunotherapy on urothelial cancer patients in data set IMvigor210at all classification thresholds (AUC=0.604). ns, not significant; *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001.

Figure 8

Pan-cancer analysis of the correlation between risk score and immune cell infiltration as well as patients’ survival. (A) Pan-cancer analysis of risk score in 30 non-hematologic tumors. (B) Pan-cancer analysis of risk score indicated patients with high risk (above median) had higher death rate (shown in pink). (C–E) Pan-cancer analysis showed that patients with lower risk score have significantly longer overall survival (C), longer disease-specific survival (D) and longer disease-free interval (E) than those with higher risk score. (F) Correlation analysis of risk score and immune cell infiltration level in 30 non-hematologic tumors by TIMER.