Antigen presentation and tumor immunogenicity in cancer immunotherapy response prediction

1Immunotherapy, represented by immune checkpoint inhibitors (ICI), including anti-PD-1 antibodies, anti-PD-L1 antibodies, anti CTLA-4 antibodies or their combinations, is transforming the treatment of cancer. Compared to conventional therapies, ICI can induce significantly improved clinical responses in patients with various types of late-stage metastatic cancers. However, the majority of unselected patients will not respond to ICI. Most tumor types show response rates below 40% to PD-1 inhibition, and the response rates of each tumor type are reported to be correlated with the tumor mutational burden (TMB) of that tumor type. Multiple factors are reported to affect ICI effectiveness, including: PD-L1 expression, TMB, DNA mismatch repair deficiency, the degree of cytotoxic T cell infiltration, mutational signature, antigen presentation defects, interferon signaling, tumor aneuploidy and T-cell signatures. These biomarkers have various rates of accuracy and utility, and the identification of a robust ICI-response biomarker is still a critical challenge in the field.

2ICI help a patient’s immune system to recognize and attack cancer cells. The immunogenicity of cancer cells is the fundamental determinant of ICI response. Theoretically, tumors of very low or no immunogenicity will not respond to therapeutic strategies that enhance the immune response. Hence, ICI can only be used to treat tumors that have sufficient immunogenicity. Furthermore, enhancing tumor immunogenicity can potentially transform an immunotherapy-non-responsive tumor into an immunotherapy-responsive tumor.

3The actual immunogenicity of a tumor is not easy to measure. In theory, tumor immunogenicity is determined by the tumor cell itself, and is also influenced by factors related to the tumor microenvironment, such as the functioning of professional antigen-presenting cells like dendritic cells (DCs). Fundamental determinants of tumor immunogenicity include tumor antigenicity, and antigen processing and presenting efficiency.

4Antigen presentation defects have already been shown to contribute to ICI-response failure. To measure antigen processing and presenting efficiency systematically, we applied a gene set variation analysis (GSVA) method to generate an antigen processing and presenting machinery (APM) score (APS), which was calculated from the mRNA expression status of APM genes. Tumor immunogenicity score (TIGS) was then calculated by combining the APM score and the TMB. The antigen-presentation gene expression signature and tumor immunogenicity landscape of 32 cancer types from The Cancer Genome Atlas (TCGA) project are provided.

5TIGS exhibits improved performance in both pan-cancer ICI objective response rate (ORR) correlation and accuracy of ICI clinical response prediction when compared with TMB. Our results suggest that TIGS represents a novel and effective tumor-inherent biomarker for the prediction of immunotherapy response.

6Cell surface presentation of peptides by major histocompatibility complex (MHC) class I molecules is critical to CD8+ T-cell mediated adaptive immune responses, including those against tumors. The generation and loading of peptides onto MHC class I molecules require the functioning of the APM. Several steps are involved in this process, including: 1) peptide generation and trimming in the proteasome; 2) peptide transport; 3) assembly of the MHC class loading complex in the endoplasmic reticulum (ER); and 4) antigen presentation on cell surface.

7The efficiency of antigen processing and presentation is one determinant of tumor immunogenicity. Here, we used the mRNA expression status of genes involved in the APM process as an indicator of the efficiency of these antigen-processing and -presenting steps. A GSVA approach was applied to measure the overall expression enrichment of APM genes. On the basis of a review paper about APM, the following genes were selected for quantification: PSMB5, PSMB6, PSMB7, PSMB8, PSMB9, PSMB10, TAP1, TAP2, ERAP1, ERAP2, CANX, CALR, PDIA3, TAPBP, B2M, HLA-A, HLA-B and HLA-C (Figure 1—source data 1). GSVA calculates the per sample overexpression level of a particular gene list by comparing the ranks of the genes in that list with those of all other genes. The resulting GSVA enrichment score is defined as the APS.

8To explore the pan-cancer distribution pattern of APS, we analyzed about 10,000 tumors of 32 cancer types from TCGA (Figure 1). The boxplot in Figure 1A shows large variance in APS across TCGA cancer types, which uncovers significant distinction in antigen-processing and -presenting efficiency among different cancer types. This analysis is similar to a previous study of seven APM genes whose expression signature is highly correlated with the APS quantified in this study (Figure 1—figure supplement 1). Patient Harmonic Best Rank (PHBR) I and II scores have recently been proposed to quantify a patient’s antigen presentation ability on the basis of the genotypes of their MHC class I or class II genes, respectively. However, no significant correlations can be observed between APS and PHBR scores (Figure 1—figure supplement 1), probably because these two methods capture different information about antigen presentation: PHBR are based on MHC genotype information, whereas APS are based on information about the expression of antigen-presentation genes. Univariate Cox regression analyses suggest that APS is associated with cancer patients' survival, and some are statistically significant (Figure 1B). Meta-analysis with pan-cancer hazard ratio values suggests that APS do not associate with prognosis (Figure 1B).

9To identify the specific gene signatures that determine patients’ APS status, we initially ran differential gene expression analysis for each TCGA cancer type on the basis of APS status. Patients with APS above the median were defined as ‘APS-High’, patients with APS below the median were defined as ‘APS-Low’. Differential expression genes (p-value < 0.01, FDR < 0.05) were ranked by logFC from high to low and then selected for gene set enrichment analysis (GSEA) with gene sets from MSigDB. In results from hallmark gene sets, several gene signatures (especially interferon alpha/gamma response) were found to be enriched in most TCGA cancer types with high APS, suggesting that high APS is strongly associated with the interferon alpha/gamma signaling pathway (Figure 2A). GSEA using Reactome gene sets further validated this result (Figure 2—figure supplement 1). Interestingly, interferon gamma was reported to regulate APM gene expression, which is consistent with this observation.

10Immune infiltration score (IIS) was calculated with GSVA using a list of marker genes for immune cell types and has been validated by the CIBERSORT method (Figure 2—source data 1). TIMER is another method that can accurately resolve the relative fractions of diverse cell types on the basis of gene expression profiles from complex tissues. To further validate the calculated IIS, we performed TIMER analysis and found that the TIMER results were highly correlated with the calculated IIS (Figure 2—figure supplement 2). Significant associations between APS and IIS at both the level of cancer types and the level of individual patients were observed (Figure 2B and C). The gene list for APS calculation did not overlap with the gene list for IIS calculation.

11Pan-cancer distribution of TMB was also analyzed with the TCGA dataset (Figure 2—figure supplement 3). Different cancer types show different prognosis in relation to high TMB (Figure 2—figure supplement 3). Meta-analysis including all TCGA cancer types suggests that patients with high TMB tend to have poor prognosis (Figure 2—figure supplement 3). TMB reflects tumor antigenicity and predicted improved survival after immunotherapy. However, in cancer patients not treated with immunotherapy, high TMB tends to be associated with poor prognosis, probably because tumors accumulate mutations during progression as a result of genome instability, and consequently, high TMB is usually associated with late-stage cancer.

12The immune cell subsets were assessed with both IIS and CIBERSORT methods, and the associations between immune cell subsets with APS were analyzed further (Figure 2—figure supplement 4). Several types of immune cells, including cytotoxic cells, show strong correlation with APS values (Figure 2—figure supplement 4). TMB and IIS show relatively weak intercorrelation (Figure 2D and E). The significant correlation between APS and IIS could be due to the following reasons: first, the immune response coordinated by interferon signaling could regulate both APS and IIS; and second, the immunogenicity contributed by APS could stimulate immune response.

13Tumor immunogenicity is determined by two factors: the antigenicity of tumor cells and the processing and presentation of tumor antigens. These two factors are independent, and are both required for tumor immunogenicity determination. Theoretically, tumor immunogenicity score (TIGS) can be represented as [“Tumor antigenicity”] x [“Antigen processing and presenting status”].

14Non-synonymous tumor mutation and, consequently, the production of neoantigens can elicit immune response. Pan-cancer TMB distribution was analyzed, and log-based TMB values were found to show a Gaussian distribution (Figure 4—figure supplement 1). In addition, a previous study had already indicated that log(TMB) shows linear correlation with pan-cancer immunotherapy ORR. Thus, we used log(TMB) as a simple representation of ‘Tumor antigenicity’. APS calculated on the basis of GSVA range from −1 to 1. To multiply with tumor antigenicity, we used normalized APS values, which range from 0 to 1, as a representation of ‘Antigen processing and presenting status’.

15We calculated tumor immunogenicity score (TIGS) by using the following formula:

16TIGS were calculated for TCGA samples for which both TMB and RNA-seq gene expression data are available (32 cancer types, 8413 samples) (Figure 3A). Cancer types with high TIGS include: skin cutaneous melanoma (SKCM), diffuse large B-cell lymphoma (DLBC), colon adenocarcinoma (COAD), head and neck squamous cell carcinoma (HNSC) (Figure 3A). Univariate Cox regression analysis suggests that TIGS is associated with cancer patients' survival, and this association is statistically significant for some cancer types (Figure 3B). Meta-analysis involving all TCGA cancer types suggested that high TIGS tends to be associated with a poor prognosis in patients not treated with immunotherapy (Figure 3B), which may be due to a mechanism that is the same as that which leads to high TMB.

17Previous studies have shown that TMB can predict pan-cancer ICI ORR. Here, we evaluated and compared the performance of APS, TIGS with TMB in pan-cancer ICI ORR correlation. The ORR for anti–PD-1 or anti–PD-L1 therapy were plotted against the corresponding median APS, TIGS, TMB across multiple cancer types. In an extensive literature search, we identified 25 tumor types or subtypes for which ORR data are available. For each tumor type, we pooled the response data from the largest published studies that evaluated ORR. We included only studies of anti–PD-1 or anti–PD-L1 monotherapy that enrolled at least 10 patients who were not selected for PD-L1 tumor expression. (Identified individual studies and references are available in Figure 4—source data 1 and Figure 4—source data 2.)

18To calculate TIGS, two different approaches can be applied. In the first approach, the APS and TMB information are obtained from different studies. This approach can include a greater number of different cancer datasets. In a second approach, all APS and TMB information is obtained from the same TCGA datasets, and in this case, fewer cancer types are available for investigation. When using the first approach, in order to calculate TIGS, the median TMB for each tumor type was obtained from a validated comprehensive genomic profiling assay that was performed and provided by Foundation Medicine. The APS information for 23 tumor types was calculated on the basis of TCGA datasets, whereas the APS for Merkel cell carcinoma, cutaneous squamous cell carcinoma and small-cell lung cancer were calculated on the basis of GEO microarray datasets. Significant correlations between APS, TMB, TIGS and the ORR were observed (Figure 4). The correlation coefficients between APS and ORR and between TMB and ORR were 0.42 (p=0.038) and 0.71 (p=6.8e-5), respectively (Figure 4), suggesting that 18% and 50% of the difference in the ORR across cancer types could be explained by APS and TMB, respectively. The correlation coefficient between TIGS and ORR is 0.78 (p=5.4e-6) (Figure 4C), indicating that 60% of the difference in ORR could be explained by TIGS. These pan-cancer ORR analyses imply that TIGS performs better than TMB or APS in correlations with immunotherapy ORR. When using the second approach for TIGS calculation, TIGS still outperformed both TMB and APS in pan-cancer ORR correlation (Figure 4—figure supplement 1).

19Compared with TMB and APS, TIGS showed improved correlation with immunotherapy ORR in various types of cancer. Here, we further evaluate the performance of TIGS in predicting ICI clinical response for individual cancer patients. Recently, several prediction biomarkers for immunotherapy response that are based on gene-expression profiling have been reported. reported an IFN-γ-related mRNA expression signature that predicts clinical response to PD-1 blockade. recently described two distinct interferon-related gene expression signatures: ISG.RS, which is associated with resistance to ICI, and by contrast, IFNG.GS, which is derived from an IFNG hallmark geneset and associated with response to ICI. reported a T-cell dysfunction and exclusion gene expression signature (named ‘TIDE’ in the original paper) as a biomarker for cancer immunotherapy response. TIDE outperforms known immunotherapy biomarkers — TMB, PD-L1 expression, and interferon gamma gene expression signature — in predicting the response to immunotherapy in melanoma and lung cancer. The predictive power of TIGS in ICI clinical response was evaluated and compared with those of TMB and biomarkers based on gene expression profiling using ICI datasets, which contain both TMB and transcriptome data for individual patients. In total, two melanoma datasets and one urothelial cancer dataset were available for this analysis.

20To evaluate performance in predicting clinical response to ICI, we used the receiver operating characteristic (ROC) curve to measure the true-positive rates against the false-positive rates at various thresholds of TMB, TIDE or TIGS values (Figure 5A–C). When compared to the widely used ICI-response biomarker TMB, TIGS consistently achieved better performance in all three ICI datasets (Figure 5A–C). The predictive power of TIGS was comparable to that of TIDE in the two melanoma datasets. However, TIDE failed to predict response to immunotherapy in urothelial cancer, so TIGS showed better performance in the urothelial cancer dataset (Figure 5C). TIGS also outperforms other immunotherapy biomarkers that are based on gene expression profiling, including IIS, IFNG, ISG.RS, IFNG.GS and CD8, in all three datasets (Figure 5D–F and Figure 5—figure supplement 1). The list of genes used to calculate IFNG, ISG.RS, IFNG.GS and CD8 signatures are available in Figure 5—source data 1. Interestingly, APS itself also shows improved or similar prediction power when compared to other gene-expression-profiling-based biomarkers (Figure 5D–F and Figure 5—figure supplement 2). The expression profiles of randomly selected genes (named ‘APSr’ in Figure 5D–F), which were used as a negative control, failed to predict immunotherapy response in all three datasets.

21In all three available datasets, Kaplan–Meier overall survival curves were further compared in patients with high vs low TIDE, TMB or TIGS level (Figure 5G–O). Patients with TIGS above the median were defined as ‘TIGS-High’ while the remaining patients were defined as ‘TIGS-Low’. ‘TMB-High’, ‘TMB-Low’, ‘TIGS-High’ and ‘TIGS-Low’ were similarly defined. Comparison of survival curves showed better survival for TMB-High patients than for TMB-Low patients in all three ICI datasets, even though the difference did not reach significance in any of the three datasets, probably because of the limited sample size (Figure 5G–I). As defined in the original paper, TIDE-Low indicates low tumor immune dysfunction and low immune escape, and consequently high immunotherapy response. In the melanoma dataset, significantly improved survival was observed in TIDE-Low patients when compared to TIDE-High patients (Figure 5M). In the urothelial cancer dataset, TIDE-Low patients did not have the expected immunotherapy response (Figure 5O). However, TIGS-High patients showed significantly better survival curves than TIGS-Low patients in all three ICI datasets (Figure 5J–L). These analyses suggest that in all three available datasets, TIGS outperforms TMB and other biomarkers that are based on gene-expression profiling (TIDE, IFNG etc.) in accurately predicting clinical response to immunotherapy and in pan-cancer applicability.

22Immunogenicity is an important inherent feature of tumor cells. This feature is determined by the tumor cell itself, and is also influenced by the tumor microenvironment. Two key determinants of tumor immunogenicity are tumor antigenicity and the ability to present such antigenicity. Here, we proposed an initial method to measure the immunogenicity of a tumor. This measured tumor immunogenicity score (TIGS) shows consistently improved correlations with immunotherapy ORR in various types of cancer when compared to TMB. TIGS also shows improved performance in ICI clinical response prediction when compared with TMB and other biomarkers that are based on gene expression profiling (TIDE, interferon gamma signature and so on) in both prediction accuracy and pan-cancer applicability. Furthermore, our tumor-immunogenicity-based biomarker could guide the treatment to transform some ICI-non-responsive tumors into ICI-responsive tumors. Stimulating the APM pathway could enhance tumor immunogenicity, and possibly ICI responsiveness.

23Our study demonstrates that TIGS is an effective biomarker for ICI-response prediction. TIGS capture two key aspects of tumor immunogenicity, antigen presentation and tumor antigenicity, which could be the reason for its improved performance in ICI-response prediction when compared to known biomarkers. Furthermore, our formula for TIGS calculation can point to a new way to transform some ICI-non-responsive tumors into responsive tumors by enhancing the tumor immunogenicity. One approach is to enhance the efficiency of antigen presentation. Our GSEA indicates that interferon signaling is the top gene signature associated with APS-High, and interferon signaling has been reported to influence APM gene expression. We may enhance antigen presentation by stimulating interferon signaling in patients who are initially not responsive to ICI, especially in cancer types that have low APS, such as prostate cancer and breast cancer.

24Our study identified several cancer types in which antigen presentation status makes a significant contribution in ICI response. Breast cancer and prostate cancer have usual TMB but fairly low ICI-response rates, probably because of low APS; renal clear cell carcinoma has good ICI response rate, possibly as a result of high APS. Furthermore, our linear correlation formula — ORR = 21.4 × TIGS – 2.7 (this formula is based on the data in Figure 4C) — can be used to make hypotheses with respect to the ORR in tumor types for which anti–PD-1 therapy has not been explored. For example, we anticipate a clinically meaningful ORR of 12.3% (95% confidence interval [CI], 8.8% to 15.8%) for uterine corpus endometrial carcinoma (UCEC) on the basis of a median TIGS of 0.7.

25This study reports the first quantification of tumor immunogenicity. Several situations need to be considered for future improvement of this quantification. First, other factors including tumor germline antigen, copy number variation status, tumor purity and intra-tumor heterogeneity should also be considered to enable more accurate measurement of the antigenicity of tumor cells. Second, for quantifying antigen presentation efficiency, APM protein expression and function assessment will be more accurate than APM mRNA expression measurement. Third, other factors that influence TIGS should also be considered, including the function of professional antigen presentation cells (dendritic cells for example) in the immune microenvironment.

26This manuscript primarily focused on the cytosolic or endogenous neoantigen presentation pathway mediated by MHC class I. This does not mean that the potential neoantigen presentation by MHC class II is not important, and further studies are needed to improve the methods for the quantification of antigen presentation in cancer patients. In addition, a sex difference in the predictive power of TMB has been reported recently in lung cancer. To explore the potential sex difference in TIGS’s predictive power, we need larger datasets with more patients.

27TIGS is an extension and enhancement of the immunotherapy biomarker TMB. TIGS is tumor cell-based, and is distinct from the recent immunotherapy biomarkers immunophenoscore or T-cell dysfunction and exclusion signature. Both of these ICI biomarkers are based on tumor immune microenvironment. As a tumor inherent biomarker, TIGS can not only be used for predicting immunotherapy response, but also point ways to manipulate the immunogenicity of tumors, and consequently the response to immunotherapy.

28The pancan normalized gene-level RNA-Seq data and clinical information for 33 TCGA cohorts were downloaded from UCSC Xena (https://xenabrowser.net/) with R package UCSCXenaTools. Samples with ‘pathologic stage’ 0 or X were filtered out and only ‘sample type’ is ‘Primary Tumor’ (32 cancer types, N = 9109) were saved for further analysis. Pre-compiled, curated somatic mutations (MC3 version) for TCGA cohorts were downloaded by the R package TCGAmutations. Microarray gene expression datasets for Merkel cell carcinoma, cutaneous squamous carcinoma and small cell lung cancer were downloaded from the GEO database via R package GEOquery. Specifically, GSE39612, GSE22396, GSE36150, GSE50451, GSE99316 were identified and downloaded.

29APM gene expression status and infiltration levels for immune cell types were quantified using the GSVA method implemented in the R package GSVA. RNA-Seq or microarray datasets were provided as input and output is a near-Gaussian list of decimals that can be used in visualization or downstream statistical analysis. Lists of genes for quantifying immune cell types were as previously described. Gene lists for APM score and quantification of immune cell type are provided in Figure 1—source data 1 and Figure 2—source data 1.

30The immune infiltration score (IIS) for a sample was defined as the mean of standardized values for macrophages, DC subsets (total, plasmacytoid, immature, activated), B cells, cytotoxic cells, eosinophils, mast cells, neutrophils, NK cell subsets (total, CD56 bright, CD56 dim), and all T-cell subsets (CD8 T, T helper, T central and effector memory, Th1, Th2, Th17, and Treg cells). In vitro validation with multiplex immunofluorescence, in silico validation using simulated mixing proportions and comparison between CIBERSORT and IIS have been described previously. TIMER is another method that can accurately resolve the relative fractions of diverse cell types on the basis of gene expression profiles from complex tissues. To further validate the calculated IIS, we performed TIMER analysis and found that the result of TIMER was highly correlated with the calculated IIS (Figure 2—figure supplement 1).

31Original APM scores (APS) from GSVA are in the range of −1 to 1. To calculate TIGS, original APM score from GSVA implementation was rescaled by the minimal and maximal APM score from TCGA Pan-cancer analysis. The formula iswhere is the minimal APM score among TCGA pan-cancer samples; and is the maximal APM score among TCGA pan-cancer samples. The normalized APM scores are in the range of 0 to 1. The normalized APS is set to 0 if a loss of function mutation exists in the B2M gene.

32TMB was defined as the number of non-synonymous alterations per megabase (Mb) of genome examined. As reported previously, we used 38 Mb as the estimate of the exome size. For studies reporting mutation number from whole exome sequencing, the normalized TMB = (whole exome non-synonymous mutations)/(38 Mb).

33We calculated TIGS as following:

34The natural logarithm was used here. Notably, some tumors have a TMB level below one mutation/Mb, so to avoid a negative number in quantifying ‘tumor antigenicity’, we added a pseudo count of one to normalized TMB. So the TIGS formula is:or

35The dataset search strategy for assessment of cancer immunotherapy ORR) assessment has been described previously. We searched MEDLINE (from January 1, 2012 to September 1, 2018), as well as abstracts in the American Society of Clinical Oncology (ASCO), the European Society for Medical Oncology (ESMO), and the American Association for Cancer Research (AACR), to identify clinical studies for anti-PD1 or anti-PDL1 therapy in various tumor types or subtypes. We searched for clinical trials using the following keywords: nivolumab, BMS-936558, pembrolizumab, MK-3475, atezolizumab, MPDL3280A, durvalumab, MEDI4736, avelumab, MSB0010718C, BMS-936559, cemiplimab, and REGN2810. We excluded studies that enrolled fewer than 10 participants, studies that investigated anti-PD-(L)one therapies only in combination with other agents, and studies that selected patients on the basis of PD-L1 expression or other immune-related biomarkers. Of the remaining studies, only the largest published study for each anti-PD-(L)one therapy was included in the final assessment of pooled ORR for each tumor type or subtype. The final identified individual studies are summarized and presented in Figure 4—source data 1. The TMB information for major solid tumor types or subtypes has been described previously. The APS of most tumor types or subtypes are based on TCGA RNA-seq data, except those for Merkel cell carcinoma, cutaneous squamous carcinoma and small cell lung cancer, which do not have available TCGA RNA-seq data. For these cancer types, the GEO datasets GSE39612, GSE22396, GSE36150, GSE50451, GSE99316 were used to generate APS. In total, 28 cancer types have both TMB and ORR values, and 25 of them also have transcriptome data that can be used for calculating APS. Therefore, TIGS were calculated for these 25 cancer types which have both TMB and APS information available (Figure 4—source data 2). Linear regression models were constructed to correlate ORR with APS, TMB and TIGS for each of the cancer types or subtypes.

36To evaluate the power of TIGS to predict clinical response to ICIs, we searched PubMed for ICI clinical studies for which TMB and gene transcriptome information was available for individual patients. In total, three datasets were identified after this search. The dataset was downloaded from the supplementary files of reference. This dataset related to CTLA-4 blockade in metastatic melanoma, and defined ‘clinical benefit’ using a composite end point of complete response or partial response to CTLA-4 blockade as assessed by RECIST criteria or stable disease by RECIST criteria with overall survival greater than 1 year, ‘no clinical benefit’ was defined as progressive disease by RECIST criteria or stable disease with overall survival less than 1 year. The dataset was downloaded from the supplementary files of reference. This dataset related to anti-PD-1therapy in metastatic melanoma: responding tumors were derived from patients who have complete or partial responses or stable disease in response to anti-PD-1 therapy; non-responding tumors were derived from patients who had progressive disease. The dataset was downloaded from https://github.com/hammerlab/multi-omic-urothelial-anti-pdl1. This dataset related to PD-L1 blockade in urothelial cancer: durable clinical benefit was defined as progression-free survival >6 months. RNA-Seq data were used to calculate the APS for each patient. Only patients for whom both APS and TMB value were available were used to calculate the TIGS. The median of TMB or TIGS was used as the threshold to separate the TMB-High and TMB-Low groups or the TIGS-High and TIGS-Low group in Kaplan-Meier overall survival curve analysis.

37The immunotherapy clinical response prediction performance of TIGS and APS have been compared with those of the following biomarkers: TMB, TIDE, IFNG, IFNG.GS, ISG.RS, PDL1, IIS, and CD8. The TIDE score was calculated using online software that is available on the website http://tide.dfci.harvard.edu. We followed the instructions on the website to generate input data for TIDE score calculation and exported the results to CSV files. The TIDE scores in the result files were used to predict response. The calculation of scores for the gene-expression-profiling-based biomarkers (i.e. IFNG, CD8, and PDL1) has been described by. The average expression values among all members defined by the original publications were used to quantify each biomarker. The interferon gamma gene expression signature (IFNG) used genes IFNG, STAT1, IDO1, CXCL10, CXCL9, and HLA-DRA. The calculation of IFNG.GS and ISG.RS scores were previously described in. CD8 used genes CD8A and CD8B. PDL1 used gene CD274. As a negative control, we performed GSVA with 18 randomly selected genes, and the resulting score was named ‘APSr’ here. This GSVA with random genes was repeated for 100 times, and APSr were used to predict immunotherapy response. The average AUC of these 100 APSr is shown.

38Univariate cox analysis was performed by R package survival. P values were adjusted using the FDR method, and FDR < 0.1 is considered statistically significant. Hazard ratios and their 95% confidence intervals for TCGA cancer types were collected and used for meta-analysis with the random effect model in the R package metafor. The receiver operator characteristic (ROC) curve was generated by plotting the rate of response at various threshold settings of TMB, TIDE or TIGS within the R package pROC. The area under the curve (AUC) was reported for each analysis. On the basis of the median of TMB, TIDE or TIGS, we separated patients into High and Low group in the survival analysis. Keplan-Meier curves of overall survival were thus plotted with log-rank test p-value in the R package ggpubr. For GSEA enrichment analysis, we compared samples that had APS above the median with those that had APS below the median across TCGA tumor types using the limma package. Genes with p-value < 0.01 and FDR < 0.05 were ranked by logFC from top to bottom and then inputted into the GSEA function of the R package clusterProfiler with custom gene sets downloaded from Molecular Signature Database v6.2. Normalized enrichment score (NES) was used to rank the differentially enriched gene sets. Correlation analysis was performed using the spearman method. All reported p-values are two-tailed, and for all analyses, p<=0.05 is considered statistically significant, unless otherwise specified. Statistical analyses were performed using R (version 3.6.0).

39All of the code and data used to generate the figures are freely available at https://github.com/XSLiuLab/tumor-immunogenicity-score (; copy archived at https://github.com/elifesciences-publications/tumor-immunogenicity-score). Analyses can be read online at https://xsliulab.github.io/tumor-immunogenicity-score/. Source data files have been provided for Figures 1, 2, 4 and 5.