Construction of immune scores to predict patient prognosis and response of anti-PD1 treatment in stage III-IV melanoma

Background: Tumor-inltrating immune cells were demonstrated to be associated with patient survival and responses of targeted therapy. However, in melanoma, there is no reliable and individualized prognostic signatures based on comprehensive evaluation of the immune prole inferred from bulk tumor transcriptomes. In this study, we aimed to develop immunoscores associated with prognosis and responses of anti-PD1 targeted therapy in stage III-IV melanoma. Methods: The immunoscore and immunoscore-based prognostic nomogram were constructed based on the melanoma cohort from the Cancer Genome Atlas (TCGA), and validated in the population from the Gene Expression Omnibus (GEO). Besides, in another cohort obtained from the GEO database, we developed an immunoscore for predicting responses of anti-PD1 therapy. Twenty-two types of immune cell fraction were estimated using CIBERSORT. The least absolute shrinkage and selection operator (Lasso) regression model was utilized to develop individualized immunoscores. Results: With the Lasso regression, an immunoscore was constructed consisting of nine types of immune cell subtypes. In both of the training (192 cases) and validation (227 cases) cohorts, signicant difference was observed between immunoscore-low and immunoscore-high groups in overall survival (OS). Multivariable analysis demonstrated that the immunoscore was an independent prognostic factor (P < 0.001) for OS. The prognostic value of the immunoscore was also conrmed by the ROC curves. Nomogram integrating immunotype and other clinical characteristics also showed good discrimination, calibration and usability in both of the training and validation cohorts. Finally, in another GEO cohort (218 cases), an immunoscore was constructed based on nine immune cell types for predicting anti-PD1 therapy response. Conclusion: The proposed immunoscores represent promising models for estimating overall survival and anti-PD1 treatment response in patients with stage III-IV melanoma.


Introduction
Cutaneous malignant melanoma (abbreviated as melanoma below) is a type of highly aggressive cancer. 1 Melanoma only accounts for nearly 2% of skin cancers, but causes the most deaths of skin cancers. 2 Worldwide, cutaneous melanoma leads to 55,500 deaths annually (0.7% of all cancer-related deaths). 1 The most commonly used system for melanoma classi cation is the American Joint Committee on Cancer (AJCC) classi cation of malignant tumors (TNM). 3 Targeted therapy and immunotherapy gradually become the new standards of care in selected patients with stage III-IV melanoma. 4 Previous studies have established phenotype-speci c gene signatures to estimate patient prognosis and predict the response to targeted therapy or immunotherapy. 5 However, up to date, these classi cations failed to be widely used in clinic practice.
Melanoma has signi cant tumor and tumor microenvironment heterogeneity. 6 Given the obvious intratumor phenotype heterogeneity and phenotype plasticity in melanoma progression and therapy response, 7-9 the evaluation and prediction of prognosis and treatment response for melanoma is complicated and challenging. 10 The tumor immune response is found to be tightly associated with tumor initiation, progression, and metastasis. 11,12 Owing to the complex interaction between the immune cell in ltration and tumor cell, there were studies illustrating that special subset of T cells, such as CD8 + TCF7 + T cells, were associated with the response to anti-PD1 therapy in melanoma; 13 and regulatory T cells (Treg) in melanoma is predictive of patient survival. 14 In this study, based on the tumor-in ltrating immune cell, we aimed to explore the association between immune scores and patient prognosis/anti-PD1 therapy response. Previous studies or guidelines usually evaluated the immune cell in ltration based on immunohistochemistry staining, and then established the immune score by integrating the immune cell numeration in the tumor core and invasive margin. 12,15,16 However, it is di cult to standardize this process and the results sometimes could not be repeated. In the last decade, with the improvement of high-throughput detection and screening method, several studies have developed various types of algorithms to comprehensively illustrate the immune pro le in tumors, such as CIBERSORT and ESTIMATE. 17,18 In the present study, based on data from different databases, we constructed immune scores to achieve accurate prediction of patient prognosis and anti-PD1 therapy responses.

Patients And Methods
Data mining of TCGA cohort and GEO database The analyzed TCGA data of skin melanoma was downloaded from The Cancer Genome Atlas (TCGA; https://portal.gdc.cancer.gov). Clinical and mRNA sequencing data of patients with stage III-IV melanoma were used in the current study. The dataset used included 192 melanoma patients and the clinicopathologic data collected included age, race, sex, T-stage, lymph node status, AJCC-TNM stage (7th ), ulceration, tumor location, Clark classi cation, Breslow depth, survival status and survival months. Data of TCGA was utilized to construct the immune score for predicting patient survival.
We performed a search in the Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo/) for eligible datasets meeting the following criteria: expression pro ling by array or high throughput sequencing; more than 20 cases in the dataset; not cell lines or other irrelevant samples (such as mouse tissues); samples with available clinical data (Fig. 1). Finally, ve datasets were enrolled including GSE19234, GSE22153, GSE22154, GSE54467 and GSE78220. All 227 patients with stage III-IV melanoma from these datasets (abbreviated as GEO1) were used to validate the immune score and prognostic model integrating clinicopathologic parameters and immunotype (developed based on the TCGA data).
In addition, we searched the GEO database to identify the gene expression dataset with anti-PD1 treatment information. Except for the criteria mentioned above, patients without available anti-PD1 treatment responses were excluded (Fig. 1). In total, ve datasets (GSE115821, GSE123728, GSE78220, GSE91061 and GSE93157) were retained and utilized in the nal analysis. There were 218 cases with stage III-IV melanoma were included for analysis of the association between immune cell in ltration and anti-PD1 treatment response (abbreviated as GEO2). The detailed analytic process is shown in Fig. 1. This study was reviewed and approved by the Ethics Committee in the West China Hospital. Cibersort Estimation CIBERSORT (http://cibersort.stanford.edu/) is a type of deconvolution algorithm, which could characterize the cell composition of complex tissues according to the normalized gene expression pro les. 17 In this study, three groups of patients (TCGA, GEO1 and GEO2) were rstly used for CIBERSORT fraction score calculation. Brie y, gene expression data with standard annotation was uploaded to the CIBERSORT web portal, and the LM22 gene signature (including 547 genes related to 22 immune cell subsets) with 100 permutations was utilized to run the CIBERSORT algorithm. For each case, the sum of all output estimates of 22 immune cell type fractions equaled 1, and the cell fractions of immune cell types could be used for comparison in different datasets.

Data analysis
Least absolute shrinkage and selection operator (Lasso) COX regression model (based on the "glmnet" package) was applied for covariate screening among 22 immune cell types. A 10-fold cross validation approach was utilized to con rm the optimal values of the penalty parameter λ. In the TCGA cohort, after Lasso cox analysis, all coe cients in the nal model were used to construct the immune score. Survival analyses were conducted by Kaplan-Meier method and examined by the log-rank test. In addition, univariate analysis was carried out by the Cox regression model. Correlation analyses of immune checkpoint molecules and the prognostic immune score were conducted by Pearson's correlation test.
Binomial logistic regression penalized by the Lasso method was used to obtain a set of immune cell types that could sensitively and speci cally predict anti-PD1 treatment response.
According to results of multivariate Cox analysis, we established a nomogram (by "rms" package) to achieve a better prognosis prediction. To test the discriminate ability of the prognostic model, the timedependent receiver operating characteristic (ROC) curve was depicted by the "timeROC" R package, and area under the ROC curves (AUC) was used to evaluate the discriminative ability of the prognostic model. Besides, the calibration curves were created to compare the predicted and actual probabilities. Decision curve analysis (DCA) was also conducted (by "rmda" package) to evaluate the usability of the nomogram. Gene set enrichment analysis (GSEA) analysis was performed to investigate the potential mechanisms based on the database of c2 (c2.cp.kegg.v6.1.symbols). The k-means cluster was performed by "cluster" and "factoextra" packages. All data analyses were completed by R version 3.6.3. P values less than 0.05 were considered statistically signi cant.

Patient characteristics
In this study, we included a total of 637 patients with stage III-IV cutaneous malignant melanoma, consisting 192 cases in the training cohort (TCGA), 227 patients in the validation cohort (GEO1) and 218 cases in the anti-PD1 treatment group (GEO2). The clinicopathologic information of patients in the training and validation cohort was presented in Table 1. Based on these two groups of patients, we established an immune score for predicting patient survival. Detailed clinical characteristics of cases receiving anti-PD1 treatment were given in Table S1. As shown in Fig Table S2.
The performance of immune score in the training and validation groups Firstly, the predictive accuracy of the immunoscore was evaluated by the ROC curves at the time points 1, 3 and 5 years (assessed as a continuous variable). As presented in Fig. 3A (training cohort) and 3B (validation cohort), the immunoscore displayed acceptable predictive ability, especially in the validation cohort. In addition, with the cut-off value of -0.69, we divided patients into immunoscore-low (IS-low) and immunoscore-high (IS-high) groups (immunotype). In all of the training (Fig. 3C), validation (Fig. 3D) and total (Fig. 3E) cohorts, patients in the IS-low group showed signi cantly better OS than those in the IShigh group (all P values < 0.05). The restricted cubic spline of immune score in the training dataset showed that the immunoscore presented a linear pro le (Fig. 3F).
The univariate analysis of the immunoscore in the training, validation and total cohorts showed that both of the immunoscore and immunotype were signi cantly associated with overall survival in all three cohorts (Table S3). In addition, Table S3 also showed the results of univariate analyses of clinicopathologic variables and all 22 types of immune cells. Besides, to further con rm the relationship between immunoscore (immunotype) and patient survival, we performed subgroup analyses based on available clinicopathologic features in the total cohort. The results demonstrated that the immunoscore and immunotype were associated with patient survival in most of the subgroups (Table S4). Based on the best optimal cut-off value of the immunoscore (-0.69), we divided patients into IS-low and IS-high groups.
In the training cohort, multivariable analysis revealed that the immunotype was an independent prognostic factor for overall survival in stage III-IV melanoma (HR = 3.25, 95% CI 2.11-5.00, P < 0.001) ( Table S5).

Association of immune score with clinical parameters, molecules, and biological pathways
In the training cohort, the prognostic immunoscore for stage III-IV melanoma was found to be negatively correlated with gene expression of certain immune checkpoint regulators (CD274, CTLA4, HAVCR2, LAG3, PDCD1, TIGIT, TNFRSF4 and CD47) and effector molecules of cytotoxic activity (GZMB and IFNG) (Fig. 4A-J). In addition, we performed GSEA analysis to explore the biological functions of the immunoscore. The results revealed that genes highly expressed in the IS-low group showed signi cant enrichment in multiple immune biological processes such as antigen processing and presentation, toll like receptor signaling pathway and natural killer cell mediated cytotoxicity (Fig. 4K). Finally, in the training cohort, we analyzed the correlation between immunoscore and clinical characteristics. As shown in Fig. 4L-O, we found that the female patient, the White, melanoma without ulceration and non-primary metastatic tumor tended to have lower immunoscore.

Nomogram Based On Immune Score And Clinical Features
To provide a quantitative instrument to predict the probability of OS, a nomogram integrating both of the immunotype and clinicopathologic indicators was established based on cases from the training cohort. By using multivariate analysis, we identi ed four prognosis-relevant factors including age, race, AJCC-TNM stage and immunotype (Table S5). In the nomogram, each value or subtype of these factors was assigned a score on the point scale. After calculating the total score and locating this value on the total point scale, we could determine the estimated probability of the 2-, 3-and 5-year survival probability ( Fig. 5A; Table S6). The 2-, 3-and 5-year ROC curves were plotted to re ect prognostic accuracy the nomogram (Fig. 5B). We found that, compared to TNM stage alone, our nomogram could predict patient survival more accurately (Fig. S1). In addition, the 2-, 3-and 5-year calibration curves depicted in Fig. 5C-E demonstrated that the nomogram performed well when compared to the performance of an ideal model in the training cohort. The decision curve analysis (Fig. 5F-H) also showed that using the nomogram could add bene t to the prediction of survival of melanoma patient.
We then validated the discrimination, calibration and usability of the nomogram in the GEO1 cohort. The 2-, 3-and 5-year AUC curves, calibration curves and DCA curves were shown in Fig. S2.
Establishment of an immune score associated with response of anti-PD1 treatment Finally, in the GEO2 cohort, the correlation analysis showed that the prognosis-relevant immunoscore could not accurately predict the anti-PD1 treatment responses (data not shown). After k-means clustering of the 22 types of immune subsets, patients were divided into two groups with different immune pro les (FIG. S3A-B). However, the two clusters of patients also did not show signi cant distinctions in response rate of anti-PD1 treatment in stage III-IV melanoma (FIG. S3C). The above results illustrated that neither the prognosis-relevant immunoscore nor the whole immune landscape diversity could accurately predict the response status to anti-PD1 therapy.
Based on the above results, we next explored whether some speci c immune subsets were associated with the response to anti-PD1 therapy. We used the Lasso logistic regression to identify immune subsets for predicting the response to anti-PD1 therapy (FIG. S3D-E). Finally, several types of immune cells were identi ed, and an immune score was constructed based on these immune cell subsets. The immune score predicting anti-PD1 therapy response is as follows: -2.665 × B cells naïve + 2.450 × plasma cells − 10.711 × T cells follicular helper − 5.721 × T cells regulatory + 1.022 × T cells gamma delta − 0.824 × monocytes + 1.435 × macrophages M0-10.820 × dentritic cells activated + 2.728 × mast cells resting. The AUC for the immunoscore is 0.78 (FIG. S3F).

Discussion
In this study, an immune score and an immunotype-based nomogram were constructed to estimate longterm overall survival in stage III-IV melanoma. The predictive ability of the nomogram for survival was superior to that of TNM stage, which was predominantly attributed to the effect of the immune score on patient prognosis. We validated the predictive values of the immunoscore and nomogram in an independent patient population from the GEO database. Besides, in the present study, we found that the immune cell type was distinct for immunoscores associated with patient survival and anti-PD1 treatment response.
Compared to stage I-II melanoma, there is tumor spread to the regional lymph nodes or more distant organs for stage III and IV melanoma. 19 Patients with advanced malignant melanoma (stage III-IV) have signi cantly poorer prognosis with 5-year survival rates for stage-III tumor ranging from 40%-78% and stage IV tumor ranging from 9%-21% in the previous study. 20 Accurate long-term survival prediction for these patients was crucial for clinical decision-making and treatment assignment. Traditional staging systems mainly focused on the tumors, without considering host immune response. 3 Given the crucial role of immune response in tumor control, it is time to begin to incorporate the immunoscore as a prognostic indicator and a component of tumor classi cation. 12,[21][22][23] In this study, we predicted the long-term survival based on the immunoscore calculated by the CIBERSORT algorithm, with the high-through sequencing and array data. Utilizing this computational analytical modality to public gene expression data, it was feasible to perform a precise investigation of more immune cell subtypes within a large patient cohort. The immunoscore enrolled several immune cell types that were signi cantly related to patient survival, like activated NK cells and CD8 T cells, and they are also the main force for anti-tumoral immune responses. [24][25][26] Based on the LM22 gene signatures, we could identify some special immune cell subtypes such as M0 and M1 macrophages. We found that, in melanoma, higher M1 macrophage was associated with better OS, while M0 macrophage showed an opposite effect. The function of macrophages, especially M1 and M2, has been widely reported in melanoma. 27,28 M2-polarized macrophage was found to be associated with tumor progression, while M1 macrophage showed pro-in ammatory and tumoricidal properties. 29,30 Further studies are still needed to validate the prognostic role of macrophages observed in the present study.
As a new approach for the classi cation of cancer, immunoscore was also used for prognosis prediction in other type of cancer, such as gastric cancer, colorectal cancer, liver cancer and lung cancer. [31][32][33][34] In stage III-IV melanoma, this is the rst study aiming to evaluate patient prognosis according to immunoscore in two independent cohorts. In addition, we built a nomogram to predict the overall survival for late-stage melanoma patients. Apart from immunotype and AJCC-TNM stage, prognostic factors in the nomogram included age and race. Consistent with the previous study, older age at diagnosis was signi cantly related to a poorer outcome. 35 In addition, the existed publication found that overall survival for cutaneous melanoma in the non-White was signi cantly lower (despite higher incidence of cutaneous melanoma in the White). 36 In our study, the Asians also had worse long-term survival. These results indicated that more emphasis is needed for melanoma screening and awareness in the non-White populations to improve the long-term survival.
Previous studies have reported several indicators for predicting response of anti-PD1 treatment. For example, the presence of CD8 + T cells at the tumor invasive margin was found to be a prerequisite for the therapeutic success of PD-1 blockade therapy. 37 A high tumor mutational burden was also correlated with the response to anti-PD-1 in melanoma. 38 However, these observations failed to be widely used clinically. 5 In this study, although the survival-predicting immune score was associated with multiple immune checkpoint molecules, it failed to predict anti-PD1 treatment response accurately. Based on results of the cluster analysis, we found that the two different clusters re ecting the immune in ltration landscape still could not discriminate response and non-response subgroups for anti-PD1 treatment. Finally, we developed an immune score (integrating nine types of immune cells) by Lasso Cox to estimate the anti-PD1 treatment responses, and it showed good discrimination. In conclusion, speci c tumor microenvironment immune cell activity could affect the e cacy of anti-PD1 treatment, such as regulatory T cells and γδ T cells. 39,40 The immunoscore for predicting response of PD-1-targeting therapies should be validated in larger population in further studies.
There were several limitations in this study. First, cases in stage III-IV melanoma is signi cantly heterogeneous, our immunoscore and prognostic nomogram still need to be validated in further studies; Second, samples included in the present study included both primary and metastatic tumors, which may in uence the applicability of the prognostic model; Third, this study used publicly available datasets and we could not get access to all the clinicopathologic features for each patient in the GEO database; Fourth, we utilized CIBERSORT method to assess the whole immune pro le in the transcriptome level, thus, we could not evaluate the immune landscape in the tumor core and invasive margin; Fifth, given the limited sample size, we did not validate the predictive value of the immunoscore for predicting anti-PD1 treatment, thus, it should be validated in other independent cohorts in further studies.
Taken together, this study provides an immunoscore and immunoscore-based prognostic model for overall survival prediction in stage III-IV melanoma, which may be a novel clinical tool for melanoma. In addition, we also constructed an immune score to predict the treatment outcome of patients receiving PD-1-targeting therapies in late-stage melanoma.

Synopsis
This study provides an immunoscore and immunoscore-based prognostic model for overall survival prediction in stage III-IV melanoma, which may be a novel clinical tool for melanoma. In addition, we also constructed an immune score to predict the response of patients receiving anti-PD1 therapies in latestage melanoma.