A Leukocyte Infiltration Score Defined by a Gene Signature Predicts Melanoma Patient Prognosis

Melanoma is one of the most aggressive cancers and accounts for over 160,000 new cancer cases worldwide every year (1, 2). Although recent advances in melanoma therapeutics have demonstrated some level of success, melanoma patients with metastatic lesions still exhibit poor prognosis with a five-year overall survival (OS) ranging from 45% for stage III to 18% for stage IV (1, 3). Melanoma tumors are characterized by high immunogenicity and are typically infiltrated by different types of immune cells (4, 5). Previous studies have shown that high level of immune cell infiltration was associated with favorable clinical outcomes of patients in both primary (6–9) and metastatic melanoma (10–14). Jönsson and colleagues (15) defined four molecular subtypes using gene expression data from a cohort of stage IV melanoma patients and further validated them in another cohort consisting of both stage III and IV melanoma patients. Out of them, the high-immune response subtype was characterized by the high expression of immune-related genes and was shown to be associated with good prognosis.

Immune cell infiltration is determined by many factors including the genomic landscape of tumor samples (16). Several studies have shown that the number of non-synonymous somatic mutations (i.e., somatic mutation burden) was positively correlated with the level of immune cell infiltration in multiple cancer types (4, 5, 17). Other than somatic mutations, copy-number variation (CNV) has been found to negatively correlated with immune cell infiltration across multiple cancer types, which was particularly evident in melanoma (18). In line with these observations, melanoma patients with higher somatic mutation burdens are more likely to be responsive to immunotherapy (19–21), whereas patients with higher CNV are more likely to be resistant to immunotherapy (18).

Recent clinical research supports the use of several immunotherapeutic strategies, demonstrating that immune manipulation could be used to reduce the size of melanoma lesions (1, 22, 23). The MAGE-A3 immunotherapy has been tested in melanoma, which targets the tumor-specific antigen MAGE-A3 (melanoma-associated antigen 3) to increase the recruitment of immune effector cells and antigen-presenting cells (24, 25). Despite the promising effects in a subset of patients, this therapy has only modest response rate and leads to adverse effects in some patients, such as autoimmunity (25–27). Biomarkers that predict patient response have the great potential in improving the successful application of the MAGE-A3 immunotherapy. Indeed, a gene signature consisting of 84 genes, mostly immune-related, has been shown to be predictive of patient response to the MAGE-A3 immunotherapy in metastatic melanoma (28). It might be interesting to investigate the predictive value of infiltrating immune cells.

In this study, we developed a 25-gene signature to calculate the level of immune cell infiltration in melanoma. This signature includes immune cell marker genes and genes involved in immune response. We applied this signature to several melanoma gene expression datasets and for each sample we calculated a score named leukocyte infiltration score (LIS) as a proxy of immune cell infiltration level. Our results indicate that LIS is predictive of patient survival in both metastatic (stage IV) and stage III melanomas. Its predictive values remain significant after considering critical clinical variables. In addition, LIS is correlated with treatment efficacy of the MAGE-A3 immunotherapy. In practice, LIS can be easily computed by measuring the expression levels of these 25 genes using simple platforms such as RT-PCR or a custom gene panel (29–31). This signature provides a practical assay that can be used to improve the treatment efficacy of melanoma.

RNA-seq data for SKCM samples generated by The Cancer Genome Atlas (TCGA) was used to define a gene signature indicative of the level of immune cell infiltration. Level 3 TCGA SKCM RNA-seq data were downloaded from Firehose (http://gdac.broadinstitute.org/). The dataset consists of gene expression profiles for a total of 463 tumor samples and provided RSEM normalized gene expression for 20,502 genes.

We used six additional microarray gene expression datasets to perform a series of analyses, including validation of the signature, prediction of patient survival and response to immunotherapy. All six datasets are available from the Gene Expression Omnibus (GEO) database with the following accession numbers: GSE65904, GSE54467, GSE35640, GSE22155, GSE8401, GSE19234 (13, 15, 28, 32–34). Sample sizes of these datasets are 214, 79, 56, 54, 47, and 44, respectively.

GSE65904 and GSE19234 include information on disease-specific survival (DSS) and post recurrence survival, respectively, while TCGA-SKCM, GSE54467, GSE22155 and GSE8401 datasets have OS. The MAGE-A3 immunotherapeutic treatment information was only available in the GSE35640 dataset. Stage III melanoma patients' OS information is available for GSE54467 and TCGA. A summary of the seven datasets is provided in Supplementary Table S1. GSE65904, GSE54467, GSE35640, GSE22155, GSE8401 and GSE19234 were measured on one-channel microarray platforms and downloaded as processed data matrices containing the expression levels of all probesets. Probeset level expression was then converted into gene expression level. For genes that were mapped to multiple probes, the maximum average expression of the probes was used to represent the gene expression.

Clinical information of SKCM samples examined by TCGA was retrieved from Firehose (http://gdac.broadinstitute.org). The dataset contains clinical covariates including patient age, sex, pathological stage at diagnosis, metastatic tumor sites, Breslow thickness, lymph node stage and metastatic stage. It is noted that Breslow thickness was measured at the time of diagnosis, but the tumor biopsy could be taken in the metastatic tumor from the same patient.

Genomic features of the TCGA melanoma samples were calculated on the basis of MAF files and DNA-sequencing profiles, which were downloaded from Firehose (http://gdac.broadinstitute.org). Tumor mutation burden (TMB) was represented by the sum of the number of non-silent mutations in a given TCGA melanoma sample. Copy-number data in a given TCGA melanoma sample was used to calculate the CNV burden, which measures the magnitude of the copy-number aberration in a tumor sample. The CNV burden can be mathematically represented as follows:

where C_j and |{f_j}$| represents the copy number and the size DNA segments j in a sample, m is the total number of abnormal segments in the genome and N is the size of the human genome. Each sample in the TCGA-SKCM dataset was assessed for CNV burden using this equation.

From 833 immune cells-related genes, identified in previous studies (35, 36), we selected those genes that best correlated with tumor-infiltrating immune cell levels in the melanoma microenvironment. Gene expression of each of the 833 genes in the TCGA-SKCM dataset was correlated with the ImmuneScore calculated by the ESTIMATE algorithm (37). The ImmuneScore reflected the infiltration of leukocytes, which were composed of different immune cells, in the tumor sample. Twenty-five genes with the highest correlation with the ImmuneScore were selected as the most representative genes for the ImmuneScore (Supplementary Table S2). In addition, five housekeeping genes (ACTB,GAPDH,RPLP0,GUSB,TFRC) were included for normalization purposes.

We applied the 25-gene signature to calculate a score that indicates the overall infiltration level of immune cells in a tumor sample based on the expression of these 25 genes. We denoted this score as LIS. Specifically, LIS score was calculated in the following ways. First, we calculated the relative expression of these 25 genes by normalizing their original expression values against the average expression of the five housekeeping genes. Second, we calculated the LIS score as the linear combination of these genes using the following equation:

where X_i refers to the relative expression level of gene i in the sample and β_i refers to the corresponding coefficient. Particularly, the coefficient β_i for the i-th gene was calculated on the basis of the principal component analysis (PCA; ref. 38)of the relative expression values (log₂ transformed fold change of genes again the average expression of the 5 housekeeping genes) of the 25 signature genes in the TCGA SKCM RNA-seq data. Particularly, β_i is the loading of the i-th genes in the principal component 1 of using these 25 genes for PCA.

A Cox proportional hazards regression model was used to investigate the effectiveness of patient-specific LISs in predicting patients' survival (DSS, post recurrence survival, or OS). Patient samples were stratified into two groups by using an indicator function |I( {LIS \ge t} )$|⁠, where t was a pre-specified threshold. Typically, t was set to be the median LIS in the samples. We then fitted a univariate Cox proportional hazards regression model to determine the association between the dichotomized LIS and patient survival. Multivariate Cox proportional hazards regression model was used to estimate the effect of LIS on survival outcome after controlling for age, sex, and tumor stage. To compare survival between groups, Kaplan–Meier plots were used for visualization. The difference between the survival times of different groups was compared by a log-rank test.

The Kaplan–Meier estimator was implemented in the “survival” R package (version 3.4.3). Specifically, the “coxph” function was used to construct Cox proportional hazards regression models. The “survfit” function was used to generate Kaplan–Meier survival curves. The “survdiff” function was used to statistically compare the difference between survival curves.

To define a prognostic gene signature that would capture immune cell infiltration in melanoma, we conducted a series of analyses summarized in Fig. 1. First, we defined a 25-gene signature from the TCGA-SKCM dataset and established a LIS as a measure of immune cell abundance in melanoma samples (See Materials and Methods). We then pursued four directions to examine the function of LIS in melanoma; (i) its prognostic role in metastatic melanoma, (ii) the prognostic role of LIS in stage III melanoma, (iii) the prediction of responsiveness to MAGE-A3 immunotherapy treatment, and (iv) the association between LIS and genomic aberrations, including CNV and total mutation burden.

We first examined the expression of previously reported 833 immune-related genes in the TCGA-SKCM dataset. These genes included immune-stimulatory, immune-inhibitory and immune cell type biomarker genes, in the tumor samples (Supplementary Table S2). For example, CD27, CD29, and CD40 were included in the immune-stimulatory category, whereas CTLA4, PDCD1 and IL10 were in the immune-inhibitory category. CD4, CD8 and CD19 were included in the immune cell type biomarker genes. Then, for each pair of immune-related genes, we calculated their Spearman correlation and found most of the genes were positively correlated with each other, regardless of whether they were annotated as immune-stimulatory or inhibitory (Fig. 2A). On the basis of this observation, we reasoned it would be feasible to use the expression of a small subset of immune-related genes to more practically estimate the abundance of infiltrating immune cells in tumor samples.

Out of the 833 immune-related genes, we defined a 25-gene signature to estimate immune cell infiltration levels in melanoma samples (see Materials and Methods; Supplementary Table S2). The number of genes in the signature was selected on the basis to balance between cost and accuracy. Examples of genes included in the signature are CD3, CD27, CD37, CD53, CD2, CD7, IL32, IL21 and CCL5, which are known markers for immune cell infiltration and immune response. The reduced number of genes in the signature, 25 versus 833, allows for the assessment of immune cells by less expensive low throughput methods (e.g., RT-PCR) and can be readily used in a clinical setting. On the basis of the expression of this 25-gene signature, we calculated a LIS for each tumor sample, which represents the extent of overall immune cells in the corresponding sample.

To confirm whether LIS accurately reflected the immune cell infiltration in melanoma, we first calculated the LIS for TCGA melanoma samples and assessed its correlation with the ImmuneScore. We found that LIS and the ImmuneScore are highly correlated (Spearman correlation coefficient SCC = 0.95, Supplementary Fig. S1A), even after stratification between primary and metastatic melanoma samples (SCC = 0.92, Primary; SCC = 0.94, Metastatic, Fig. 2B and C).

In a recent study, Saltz and colleagues (39) used deep learning method to estimate the infiltration of leukocyte and lymphocyte, which is an important subtype of leukocyte, in the TCGA tumor samples from the IHC staining data. We first confirmed that the ImmuneScore reflected immune cell infiltration in the melanoma samples by correlating the ImmuneScore with the leukocyte percentage across all TCGA melanoma samples and found a correlation coefficient of 0.85 (SCC = 0.85, Supplementary Fig. S1B). Moreover, to further validate that the LIS reflected the presence of the leukocyte infiltration in the melanoma, we examined the correlation between LIS and the leukocyte percentage in all TCGA melanoma samples and found a strong positive correlation (SCC = 0.80, Supplementary Fig. S1C). Moreover, we also examined the correlation between LIS and the lymphocyte percentage in all TCGA melanoma samples and found a correlation of coefficient equals to 0.49 (Supplementary Fig. S1D). In addition, the TCGA consortium has clinically quantified the lymphocyte infiltration level in all melanoma patients (40). To provide additional validation that the LIS also captured the lymphocyte infiltration, we stratified all TCGA melanoma samples into tertiles based on clinically quantified lymphocyte infiltration (Lymphocyte Score; i.e., low, medium and high). We found LISs were significantly different among the three groups and captured lymphocyte abundance (P < 2E−16, Supplementary Fig. S1E). We then separated primary and metastatic samples and found identical results; LIS captured the trend of the clinically determined Lymphocyte Score in both primary and metastatic melanoma (P = 0.02, Primary; P < 2E−16, Metastatic, Fig. 2D and E). In addition, correlating LIS directly with lymphocyte infiltration in primary and metastatic tumors also showed a significant positive association (SCC = 0.30, Primary; SCC = 0.56, Metastatic, Supplementary Fig. S1F and S1G). Thus, these observations suggested that LIS reflected the leukocyte infiltration and also captured lymphocyte infiltration in TCGA melanoma samples.

To further validate LIS as a predictor of immune infiltration, we examined the correlation of LIS with the ImmuneScores in six additional datasets with microarray-based gene expression. Across all datasets, we detected a mean correlation coefficient of 0.89 (Supplementary Fig. S1I–S1L). As an example, in the GSE65904 dataset, which had the largest sample size, we observed correlation coefficients of 0.86 and 0.92 in primary and metastatic samples, respectively (Fig. 2F and G). These findings indicated that LIS reflected immune cell infiltration in the melanoma TCGA RNA-seq and other independent microarray datasets.

Because immune cells correlate with prognosis in both primary and metastatic melanoma (6–14, 41), we next examined the prognostic role of LIS in the RNA-seq TCGA-SKCM dataset. We hypothesized that high LIS, indicating the presence of immune cells, would be associated with favorable prognosis in this melanoma cohort. Indeed, LIS was prognostic in TCGA melanoma dataset, where patients with a high LIS had better OS than patients with low LIS (HR = 0.65, P = 0.007; Fig. 3A). In addition, primary melanoma samples presented with a significantly lower LIS than metastatic melanoma samples (P = 2E−6, Fig. 3B).

We continued the characterization of LIS in TCGA by investigating the association of LIS and tumor sites in greater detail. Among primary tumors, samples were classified on the basis of their anatomic sites (i.e., trunk, extremities, and head and neck). No significant differences were observed among the groups (Fig. 3C). However, grouping metastatic melanoma samples according to metastatic sites, did yield significant differences in LIS; lymph-node metastatic samples demonstrated significantly higher LIS than regional skin or soft tissue metastatic and distal metastatic samples (P = 2E−6, Fig. 3C). In addition, the lymph-node metastatic samples or regional skin or soft tissue metastatic samples have a significantly higher LIS than primary samples (All P < 0.05, Fig. 3C).

Several clinical factors are known to be related with prognosis in melanoma. For example, the AJCC staging system is a strong predictor of melanoma prognosis (42). We therefore performed a similar association analysis to identify immune cell infiltration related clinical factors. One of the important prognostic clinical factors in AJCC staging system is the Breslow thickness, a measure of the vertical depth of the primary melanoma lesion (1). Although Breslow thickness is only measured at the primary site, we correlated LIS at metastatic sites with Breslow thickness at the primary site to see if Breslow thickness was associated immune cell infiltration in metastatic lesions. LIS was negatively correlated with Breslow thickness in both primary and metastatic melanoma samples (SCC = −0.37 and −0.24 in primary and metastatic samples, respectively; Fig. 3D and E). An average correlation of SCC = −0.32 was observed using all TCGA melanoma samples, indicating an inverse association of LIS with tumor Breslow thickness (Supplementary Fig. S2A). We also tested the association between LIS and lymph node and metastatic stages, but no statistically significant associations in either primary or metastatic samples were observed (Supplementary Fig. S2B–S2E). Taken together, we concluded that patients having high Breslow thickness tended to have low LIS in both primary and metastatic tumor samples.

Efforts by the TCGA consortium classified SKCM into multiple subtypes based on genomic and transcriptomic profiles (40). To follow up with the reported heterogeneity of this disease, we first examined the association between LIS and mutational subtypes, which included mutant BRAF, mutant RAS, mutant NF1, and Triple-WT (wild-type) groups (40). This initial assessment did not detect statistically significant associations (Supplementary Fig. S2F and S2G). We then tested whether LIS and transcriptome profiles-based subtypes (keratin, MITF-low, and immune) were correlated. Patients with high LIS were more likely to be clustered into immune-related clusters, which were indicative of good prognosis (P = 1E−6, Primary; P < 2E−16, Metastatic, Fig. 3F and G).

Proto-oncogene LCK, a known biomarker of immune cells, was found to be prognostic in melanoma (40). To examine whether LIS remains prognostic after adjusting for LCK protein expression, we stratified metastatic patients into high and low LCK protein expression groups (i.e., LCK-Hi and LCK-low). For each group, the patients were further separated into LIS-Hi and LIS-Lo groups. LIS remains prognostic in LCK high group but not prognostic in the LCK low group (Supplementary Fig. S3A and S3B).

After characterizing the relation between LIS and established prognostic factors and classifications in primary and metastatic melanoma, we then investigated the ability of LIS to specifically predict patient survival in metastatic melanomas. Consistent with findings in the mixed cohort of primary and metastatic tumors (Fig. 3A), metastatic patients with higher LIS had better prognosis (HR = 0.63, P = 0.007, Fig. 4A). However, we did not detect an association between LIS and prognosis in primary melanoma samples (Supplementary Fig. S3C). We note here that among the primary tumors, there were only 10 events, a rather small number that might have contributed to the lack of statistical significance. This may have been due to the lack of complete clinical follow up.

To determine whether LIS could consistently predict prognosis in metastatic melanoma, we extended our analyses into four publicly available metastatic melanoma datasets (GEO accession number GSE65904, GSE22155, GSE8401 and GSE19234). In all four datasets, hazard ratios (HR) were generally consistent, as expected (Fig. 4B–E). For example, high LIS indicated good prognosis in the GSE65904 dataset, which had the largest sample size (HR = 0.58, P = 0.007, Fig. 4B). To investigate whether LIS provided additional information to the established prognostic clinical factors, we applied a multivariate Cox regression model. Controlling for clinical covariates (tumor pathological stage at diagnosis, patient's age and gender), the associations remain statistically significant (all P < 0.05; Fig. 4F–H). In addition, for the datasets which didn't contain the complete clinical factors (GEO accession number GSE65904 and GSE19234), we also performed the multivariate Cox regression model and found that the associations remained statistically significant (all P < 0.05; Supplementary Fig. S2H and S2I).

Although we have examined the prognostic role of the LIS in metastatic melanoma, the majority of the patients were stage IV in which the tumor samples were collected from distant metastasis.

However, the stage III patients, whose tumor samples were mostly collected from the regional lymph nodes, had a different tumor microenvironment compared with stage IV patients. Moreover, the survival outcome of stage III melanoma patients is highly variable (42, 43), it would be useful to predict it. Thus, we investigated whether LIS could predict the stage III prognosis in melanoma. We first calculated the LIS in the GSE54467 dataset and examined its prognostic role. We found a significant protective association of LIS with survival (HR = 0.33, P = 0.009, Fig. 5A). Adjusting for clinical covariates, including pathological stage at diagnosis, age, and sex, did not substantially change the significant effect we observed (HR = 0.35, P = 0.02, Fig. 5B). We were able to replicate this finding in the TCGA stage III patient dataset (HR = 0.49, P = 0.02, Supplementary Fig. S3D).

Antigen-specific immunotherapy in the form of vaccines has been investigated in several tumor types (24, 25, 44). In melanoma, MAGE-A3 immunotherapy has been developed as a melanoma-specific immunization treatment and was tested on clinical trials for stage III patients (25). Therefore, we investigated whether LIS could predict metastatic patients' response to tumor antigen vaccine therapies. To test this, we focused on patients with tumor transcriptome profiles before receiving the treatment, and calculated LIS for each. Patients who responded to MAGE-A3 consistently showed significantly higher LIS (P = 0.01, Fig. 5C). To further verify that the LIS accurately predicts response to immunotherapy, patients were stratified into tertiles based on LIS scores (i.e., LIS-Hi, LIS-Intermediate and LIS-Lo). Compared with the mean response rate across all samples, the LIS-Hi group showed a 1.5-fold increase in response rate (P = 0.04, Fig. 5D). We further examined the prediction power of LIS by calculating its area under the curve (AUC). The AUC was 0.7 by using the LIS as the only feature (Supplementary Fig. S3E). Moreover, we applied the logistic regression model by using continuous value of LIS and binary value of LIS, which used median LIS as the cutoff value, as the predictors. As the result, the LIS was a significant predictor for MAGE-A3 immunotherapy (P = 0.02, continuous LIS; P = 0.03, binary LIS, Supplementary Fig. S3F). These results suggested that LIS is a significant predictor of immunotherapeutic response.

Immune cell infiltration was associated with genomic features such as TMB and CNV (4, 5, 17, 18, 20, 21). TMB was associated with immunogenicity in tumors, as well as response to immunotherapy (19–21, 36). Thus, we tested the association of TMB in TCGA samples with LIS. Unexpectedly, we did not observe a significant association between TMB and LIS in TCGA melanoma samples (SCC = 0.02, Supplementary Fig. S3G). In addition, no significant association between TMB and LIS was observed in either primary (Fig. 6A) or metastatic melanoma (Fig. 6B).

A recent study concluded that tumor ploidy negatively correlated with immune cell infiltration in the tumor microenvironment and disease prognosis (18). To follow-up on this finding, we investigated the correlation between the CNV and LIS in the TCGA melanoma samples. We detected a significant, inverse association between CNV and LIS, consistent with the published finding (SCC = −0.56, Supplementary Fig. S3H). Both primary tumors and metastatic tumors presented a significant negative correlation with tumor ploidy (SCC = −0.47, Primary; SCC = −0.62, Metastatic, respectively; Fig. 6C and D).

We next investigated whether individual genes with either of these two genomic alterations were correlated with LIS. In prior studies, 25 significantly highly mutated genes were identified in melanoma (40). Given this observation, we examined the differences in LIS between mutated and wildtype genes (Supplementary Tables S4 and S5). Accounting for multiple hypothesis testing, no genes' mutation status was significantly associated with LIS. We further investigated the association between CNV of those genes and LIS in primary and metastatic tumors respectively (Supplementary Tables S6 and S7). We identified four genes with CNVs associated with LIS: CDK4, KRAS, DDX3X, and CDKN2A (FDR < 0.05; Fig. 6E–H). Notably, CDK4 amplification events are associated with decreased immune cell infiltration in the primary samples. In metastatic tumor samples, amplification of oncogenic “driver,” KRAS, as well as deletion of DDX3X or CDKN2A, are associated with decreased immune cell infiltration (P = 6E−5, CDK4, P = 1E−3, KRAS, P = 7E−7, DDX3X, P = 2E−5, CDKN2A).

Tumor infiltrating immune cells exhibit characteristic clinical features and prognosis in melanoma (4–15). Although some studies suggest the prognostic role of immune cells in the melanoma tumor microenvironment (6–15), a multicohort survival analysis that systemically evaluates the prognostic role of immune cells in metastatic melanoma, especially stage III melanoma with highly variable outcomes, was lacking. Many clinical studies involving immunotherapeutic treatments aim to improve prognosis for melanoma patients, particularly stage III and metastatic patients. The association between immunotherapeutic response and the immune cell infiltration needed to be more thoroughly investigated. Despite efforts of improving prognosis of melanoma patients in the field, the precise etiology and mechanism of clinical factors and genomic features associated with immune cell infiltration in melanoma need further investigation.

In this study, we defined the LIS using the normalized expression of 25 genes to allow for clinical estimation of immune cell infiltration levels in melanoma samples. These genes were selected on the basis of their correlation with the ImmuneScore in the TCGA SKCM RNA-seq dataset. We chose a total of 25 genes, as a good balance between cost and accuracy, to reflect the immune cell infiltration. In fact, gene signatures with the top 30, 35, 40, and 50 genes achieve similar results. We then applied LIS to six independent melanoma microarray datasets generated from different platforms. Our results first indicated that LIS correlated with the ImmuneScore in these datasets (Fig. 2B and C; Fig. 2F and G; Supplementary Fig. S1G–S1K). Then, we found that LIS was able to reflect the leukocyte infiltration in the tumor samples (Supplementary Fig. S1B). Furthermore, we showed that LIS correlates with lymphocyte infiltration and a clinically measured lymphocyte score (Fig. 2D and E; Supplementary Fig. S1C; Supplementary Fig. S1E and S1F). In summary, the LIS accurately reflected the immune cell infiltration by using 25 genes. Compared with the ESTIMATE algorithm which uses 141 genes (37), our method is able to make accurate predictions using substantially fewer genes. Moreover, the traditional method for the estimation of immune cell infiltration is based on IHC staining, which is semiquantitative and labor-intensive. This 25-gene signature may be more generalizable, as it can be applied using more cost-effective, low throughput methods (e.g., RT-PCR, capture-based sequencing).

We have shown that LIS is prognostic for metastatic melanoma, including patients diagnosed with stage III disease (Fig. 3A; Supplementary Fig. S3D). We performed a multicohort prognostic analysis to examine the role of LIS in five independent microarray datasets (Fig. 4; Fig. 5A and B). The prognostic role of LIS was consistent across five independent microarray datasets and its prognostic role is independent of established clinical factors (Fig. 4; Fig. 5A and B). Furthermore, we have shown that LIS could provide additional information about the immune cell marker, LCK, in predicting the prognosis in metastatic melanoma (Supplementary Fig. S3A and S3B). Although we did not observe a prognostic role of LIS in LCK-Lo group, this may be due to the fact that LCK is associated with immune cell infiltration in melanoma. Therefore, in the immune cells poorly infiltrated samples, we lacked sufficient statistical power to predict prognosis. In addition, because our 25-gene signature also included LCK gene and the genes in our signature were highly correlated with each other, it's also possible that the rest 24 genes were not strong enough to further stratify the patients in the LCK low group into high and low risk groups. Intriguingly, we found LIS was predictive of prognosis in stage III melanoma, which is known for high variability of survival time (Fig. 5A and B; Supplementary Fig. S3D). These observations provide support for using LIS to stratify stage III melanoma patients into low- and high-risk groups.

Targeted immunotherapies have been increasingly used in the clinic. Many immunotherapies are currently being tested in clinical trials (22, 23, 25, 26). Given their known efficacy in metastatic melanoma patients and the importance of treating stage III melanoma patients, we sought to address whether our prognostic gene signature is correlated with clinical outcomes in stage III patients receiving such treatments. We examined the predictive ability of LIS with regard to the response to tumor antigen-specific vaccine, MAGE-A3, and found that LIS is associated with patients' response to MAGE-A3 therapy (Fig. 5C and D; Supplementary Fig. S3E and S3F). Patients with high LIS before receiving the therapies showed a higher response rate than those with low LIS. Our predictive analysis can be further pursued as new datasets become available in the future. Recently, the immune checkpoint blockade therapies have been used to treat melanoma and increased metastatic melanoma patients OS (20, 45, 46). In addition, many clinical trials are ongoing to test the efficacy of immune checkpoint blockade therapies in improving patient survival in melanoma. Multiple gene-expression datasets that contain gene expression associated with response to the therapy were released (20, 45). We applied our 25-gene signature in the CTLA-4 and PD-1 blockade therapy cohorts to calculate LIS for each patient and examined the correlation between LIS and patient response, however, no significant association was detected.

Our results indicate that LIS is significantly inversely associated with tumor genome ploidy whereas not significantly correlated with TMB in both primary and metastatic melanoma tumor samples (Fig. 6A–D). It is known that TMB is positively associated with neoantigen burden in melanoma, which recruits immune cells into the tumor microenvironment (36, 47). Patients with high TMB was more likely to benefit from immunotherapy and presented a prolonged survival (19–21). Meanwhile, CNV burden was reported to be negatively associated with immune cell infiltration in the tumor microenvironment and was associated with the response to the immunotherapy (18). It is still under investigation that which genomic feature contributed to the immune cell infiltration in the melanoma microenvironment. In our present analysis, we identified LIS to be significantly inversely correlated with CNV burden, but not with TMB (Fig. 6A–D). It is should be noted that CNV burden and TMB might be correlated with tumor purity due to technical reasons—more CNVs can be identified in higher purity tumor samples. Indeed, CNV burden is positively correlated with purity in both primary and metastatic TCGA melanoma with SCC > 0.4. Adjusting purity as a confounding variable will underestimate the association between CNV and LIS, because high LIS samples are typically those with low purity. Nevertheless, the association between CNV and LIS remains significant in TCGA metastatic melanoma, though much weaker, after adjusting purity.

Notably, four genes had CNV events associated with changes in immune cell infiltration. CDK4 amplification was associated with low immune cell infiltration in primary melanoma samples (Fig. 6E). This observation supports the prior finding that a CDK4 inhibitor could induce antitumor immunity (48). We also found that deletion of DDX3X or CDKN2A, known to participate in innate immunity (49), were associated with immune cell infiltration (Fig. 6G and H). On the other hand, KRAS amplification was associated with low immune cell infiltration in metastatic melanoma (Fig. 6F). Our findings provide an understanding of cancer biology in melanoma by showing which somatic mutation might influence the immune cell infiltration in the melanoma microenvironment.

In summary, we defined a 25-gene signature that captures the immune cell infiltration in melanoma tumors. We have demonstrated that LIS is associated with clinical factors, genomic features and immunotherapy response. The computational framework can be readily extended to define immune cell infiltration related signatures in other cancer types.

No potential conflicts of interest were disclosed.

Conception and design: I.P. Gorlov, C.I. Amos, C. Cheng

Development of methodology: C. Cheng

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): M. Berwick, C. Cheng

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): Y. Zhao, I.P. Gorlov, E. Hernando, C. Cheng

Writing, review, and/or revision of the manuscript: Y. Zhao, E. Schaafsma, E. Hernando, N.E. Thomas, R. Shen, M.J. Turk, M. Berwick, C.I. Amos

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): N.E. Thomas, M. Berwick, C. Cheng

Study supervision: C. Cheng

This work was supported by the National Center for Advancing Translational Sciences of the National Institutes of Health under award number KL2TR001088 and the National Institutes of Health/National Cancer Center P01 CA206980-01A1, the Center of Biomedical Research Excellence (COBRE) grant under award number GM103534, and the Dartmouth Geisel School of Medicine Start-up Fund. The authors thank Youdinghuan Chen for assistance with the manuscript review.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.