Introduction {#s1}
============

Immunotherapy, 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 ([@bib48]). Multiple factors are reported to affect ICI effectiveness, including: PD-L1 expression ([@bib14]; [@bib36]), TMB ([@bib31]; [@bib37]), DNA mismatch repair deficiency ([@bib18]), the degree of cytotoxic T cell infiltration ([@bib40]), mutational signature ([@bib26]; [@bib43]), antigen presentation defects ([@bib7]; [@bib50]), interferon signaling ([@bib1]), tumor aneuploidy ([@bib10]) and T-cell signatures ([@bib17]). 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 ([@bib28]).

ICI 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.

The 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) ([@bib25]). Fundamental determinants of tumor immunogenicity include tumor antigenicity, and antigen processing and presenting efficiency ([@bib4]).

Antigen presentation defects have already been shown to contribute to ICI-response failure ([@bib7]; [@bib50]). 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) ([@bib12]), 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.

TIGS 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.

Results {#s2}
=======

APM score definition and pan-cancer analysis {#s2-1}
--------------------------------------------

Cell 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 ([@bib19]).

The 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 ([@bib12]). On the basis of a review paper about APM ([@bib19]), 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](#fig1sdata1){ref-type="supplementary-material"}). 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.

To explore the pan-cancer distribution pattern of APS, we analyzed about 10,000 tumors of 32 cancer types from TCGA ([Figure 1](#fig1){ref-type="fig"}). The boxplot in [Figure 1A](#fig1){ref-type="fig"} 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 ([@bib35]) whose expression signature is highly correlated with the APS quantified in this study ([Figure 1---figure supplement 1](#fig1s1){ref-type="fig"}). 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 ([@bib23]; [@bib22]). However, no significant correlations can be observed between APS and PHBR scores ([Figure 1---figure supplement 1](#fig1s1){ref-type="fig"}), 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](#fig1){ref-type="fig"}). Meta-analysis with pan-cancer hazard ratio values suggests that APS do not associate with prognosis ([Figure 1B](#fig1){ref-type="fig"}).

![Analysis of antigen processing and presenting machinery (APM) score in 32 cancer types.\
(**A**) APM scores were calculated with GSVA in 32 TCGA cancer types. (**B**) Results of Cox proportional hazards regression analysis using APM score for all solid cancers. Forest plots showing log~e~ hazard ratio (95% confidence interval). Cox p-values are adjusted the with false discovery rate (FDR) method, p-values less than 0.1 are in bold. The pooled hazard ratio and p-value are generated by the random effect model. The statistical test for heterogeneity is also shown in the last column. Tumor types are ordered by median APM scores.\
10.7554/eLife.49020.005Figure 1---source data 1.APM gene list for GSVA.](elife-49020-fig1){#fig1}

APS determinants and associations in cancer {#s2-2}
-------------------------------------------

To 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 ([@bib39]). 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](#fig2){ref-type="fig"}). GSEA using Reactome gene sets further validated this result ([Figure 2---figure supplement 1](#fig2s1){ref-type="fig"}). Interestingly, interferon gamma was reported to regulate APM gene expression ([@bib2]; [@bib16]), which is consistent with this observation.

![Gene expression signatures associated with high APM score.\
(**A**) Gene sets enriched in patients with high APM score. (**B**) Significant correlation between APM score and IIS in 8949 cancer samples. (**C**) Significant correlation between APM score and IIS in different cancer types. (**D**) Correlation between TMB and IIS in 8413 cancer samples. (**E**) Correlation between TMB and IIS in different cancer types.\
10.7554/eLife.49020.011Figure 2---source data 1.Immune cell types and corresponding signature gene lists for GSVA.](elife-49020-fig2){#fig2}

Immune 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 ([@bib35]) ([Figure 2---source data 1](#fig2sdata1){ref-type="supplementary-material"}). TIMER ([@bib20]) 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 ([@bib20]) and found that the TIMER results were highly correlated with the calculated IIS ([Figure 2---figure supplement 2](#fig2s2){ref-type="fig"}). 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](#fig2){ref-type="fig"}). The gene list for APS calculation did not overlap with the gene list for IIS calculation.

Pan-cancer distribution of TMB was also analyzed with the TCGA dataset ([Figure 2---figure supplement 3](#fig2s3){ref-type="fig"}). Different cancer types show different prognosis in relation to high TMB ([Figure 2---figure supplement 3](#fig2s3){ref-type="fig"}). Meta-analysis including all TCGA cancer types suggests that patients with high TMB tend to have poor prognosis ([Figure 2---figure supplement 3](#fig2s3){ref-type="fig"}). 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.

The immune cell subsets were assessed with both IIS and CIBERSORT ([@bib27]) methods, and the associations between immune cell subsets with APS were analyzed further ([Figure 2---figure supplement 4](#fig2s4){ref-type="fig"}). Several types of immune cells, including cytotoxic cells, show strong correlation with APS values ([Figure 2---figure supplement 4](#fig2s4){ref-type="fig"}). TMB and IIS show relatively weak intercorrelation ([Figure 2D and E](#fig2){ref-type="fig"}). 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.

Tumor immunogenicity score: definition and pan-cancer profiling {#s2-3}
---------------------------------------------------------------

Tumor 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"\].

Non-synonymous tumor mutation and, consequently, the production of neoantigens can elicit immune response ([@bib34]). Pan-cancer TMB distribution was analyzed, and log-based TMB values were found to show a Gaussian distribution ([Figure 4---figure supplement 1](#fig4s1){ref-type="fig"}). In addition, a previous study had already indicated that log(TMB) shows linear correlation with pan-cancer immunotherapy ORR ([@bib48]). 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'.$${APS}_{normalized} = \frac{APS - {APS}_{pancan\_ min}}{{APS}_{pancan\_ max} - {APS}_{pancan\_ min}}$$

We calculated tumor immunogenicity score (TIGS) by using the following formula:$${TIGS = APS}_{normalized} \times log\left( TMB \right)$$

TIGS were calculated for TCGA samples for which both TMB and RNA-seq gene expression data are available (32 cancer types, 8413 samples) ([Figure 3A](#fig3){ref-type="fig"}). 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](#fig3){ref-type="fig"}). 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](#fig3){ref-type="fig"}). 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](#fig3){ref-type="fig"}), which may be due to a mechanism that is the same as that which leads to high TMB.

![Tumor immunogenicity score (TIGS) analysis in 32 cancer types.\
(**A**) Analysis of TIGS in 32 cancer types. (**B**) Results of Cox proportional hazards regression analysis using TIGS for all solid cancers. Forest plots showing log~e~ hazard ratio (95% confidence interval). Cox p-values are adjusted with the FDR method. p-values less than 0.1 are in bold. The pooled hazard ratios and the p-values were generated using the random effect model. The statistical test for heterogeneity is also shown in the last column. Tumor types are ordered by median TIGS score.](elife-49020-fig3){#fig3}

TIGS and pan-cancer ORR to PD-1 inhibition {#s2-4}
------------------------------------------

Previous studies have shown that TMB can predict pan-cancer ICI ORR ([@bib48]). 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](#fig4sdata1){ref-type="supplementary-material"} and [Figure 4---source data 2](#fig4sdata2){ref-type="supplementary-material"}.)

To 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 ([@bib5]). 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](#fig4){ref-type="fig"}). 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](#fig4){ref-type="fig"}), 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](#fig4){ref-type="fig"}), 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](#fig4s1){ref-type="fig"}).

![TIGS and predicted pan-cancer response rates to PD-1 inhibition.\
Correlation between (**A**) APS, (**B**) TMB, (**C**) TIGS and objective response rate (ORR) with anti-PD-1 or anti-PD-L1 therapy in 25 cancer types. Shown are median normalized APS (A), median number of TMB (non-synonymous mutation/MB) in log scale (B) and TIGS in 25 tumor types or subtypes among patients who received inhibitors of PD-1 or PD-L1 (C), as described in published studies for which data regarding the ORR are available. The number of patients who were evaluated for the ORR is shown for each tumor type (size of the circle), along with the number of tumor samples that were analyzed to calculate the APS, TMB or TIGS (degree of shading of the circle).\
10.7554/eLife.49020.016Figure 4---source data 1.List of citations for individual studies used in pooled analysis of objective response rate.\
10.7554/eLife.49020.017Figure 4---source data 2.Summary of pooled ORR, median TMB and median APS by tumor type or subtype.](elife-49020-fig4){#fig4}

TIGS and prediction of clinical response to ICI {#s2-5}
-----------------------------------------------

Compared 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 ([@bib1]; [@bib17]). [@bib1] reported an IFN-γ-related mRNA expression signature that predicts clinical response to PD-1 blockade. [@bib3] 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. [@bib17] 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 ([@bib17]). 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 ([@bib15]; [@bib41]) and one urothelial cancer ([@bib38]) dataset were available for this analysis.

To 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](#fig5){ref-type="fig"}). When compared to the widely used ICI-response biomarker TMB, TIGS consistently achieved better performance in all three ICI datasets ([Figure 5A--C](#fig5){ref-type="fig"}). 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](#fig5){ref-type="fig"}). 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](#fig5){ref-type="fig"} and [Figure 5---figure supplement 1](#fig5s1){ref-type="fig"}). The list of genes used to calculate IFNG, ISG.RS, IFNG.GS and CD8 signatures are available [in Figure 5---source data 1](#fig5sdata1){ref-type="supplementary-material"}. Interestingly, APS itself also shows improved or similar prediction power when compared to other gene-expression-profiling-based biomarkers ([Figure 5D--F](#fig5){ref-type="fig"} and [Figure 5---figure supplement 2](#fig5s2){ref-type="fig"}). The expression profiles of randomly selected genes (named 'APSr' in [Figure 5D--F](#fig5){ref-type="fig"}), which were used as a negative control, failed to predict immunotherapy response in all three datasets.

![TIGS predicts clinical response to ICI immunotherapy.\
(**A**) ROC curves for the performance of TMB, TIDE and TIGS in predicting anti-CTLA4 Immunotherapy response in 35 melanoma patients (dataset from [@bib41]). (**B**) ROC curves for the performance of TMB, TIDE and TIGS in predicting anti-PD-1 immunotherapy response in 27 melanoma patients (dataset from [@bib15]). (**C**) ROC curves for the performance of TMB, TIDE and TIGS in predicting anti-PD-L1 immunotherapy response in 22 urothelial cancer patients (dataset from [@bib38]). (**D--F**) AUC values of TMB, TIGS, TIDE, PDL1, immune infiltration score (IIS), interferon gamma gene expression signature (IFNG), CD8, APS and random genes as negative control for APS quantification (APSr) in the [@bib41] dataset (D), the [@bib15] dataset (E) and the [@bib38] dataset (F). The performance of a random predictor (AUC = 0.5) is represented by the dashed line. (**G,J,M**) Patients were grouped on the basis of TMB (G), TIGS (J) or TIDE (M) status. The Kaplan--Meier (KM) overall survival curves were compared between TMB-High and TMB-Low (100 patients), between TIGS-High vs TIGS-Low (35 patients) or between TIDE-High and TIDE-Low (37 patients) in the [@bib41] dataset. (**H,K,N**) Patients were grouped on the basis of TMB (**H**), TIGS (K) or TIDE (N) status. The KM overall survival curves were compared between TMB-High and TMB-Low (37 patients), between TIGS-High and TIGS-Low (26 patients) or between TIDE-High and TIDE-Low (26 patients) in the [@bib15] dataset. (**I,L,O**) Patients were grouped on the basis of TMB (**I**), TIGS (**L**) or TIDE (**O**) status. The KM overall survival curves were compared between TMB-High and TMB-Low (22 patients), TIGS-High and TIGS-Low (22 patients) or TIDE-High and TIDE-Low (25 patients) in the [@bib38] dataset.\
10.7554/eLife.49020.021Figure 5---source data 1.List of genes in the lists used for CD8, IFNG, ISG.RS and IFNG.GS signature calculation.](elife-49020-fig5){#fig5}

In 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](#fig5){ref-type="fig"}). 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](#fig5){ref-type="fig"}). As defined in the original paper ([@bib17]), TIDE-Low indicates low tumor immune dysfunction and low immune escape, and consequently high immunotherapy response. In the [@bib41] melanoma dataset, significantly improved survival was observed in TIDE-Low patients when compared to TIDE-High patients ([Figure 5M](#fig5){ref-type="fig"}). In the urothelial cancer dataset ([@bib38]), TIDE-Low patients did not have the expected immunotherapy response ([Figure 5O](#fig5){ref-type="fig"}). However, TIGS-High patients showed significantly better survival curves than TIGS-Low patients in all three ICI datasets ([Figure 5J--L](#fig5){ref-type="fig"}). 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.

Discussion {#s3}
==========

Immunogenicity 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.

Our 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 ([@bib2]; [@bib16]). 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.

Our 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](#fig4){ref-type="fig"}) --- 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.

This 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.

This 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 ([@bib45]; [@bib46]). To explore the potential sex difference in TIGS's predictive power, we need larger datasets with more patients.

TIGS is an extension and enhancement of the immunotherapy biomarker TMB. TIGS is tumor cell-based, and is distinct from the recent immunotherapy biomarkers immunophenoscore ([@bib6]) or T-cell dysfunction and exclusion signature ([@bib17]). 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.

Materials and methods {#s4}
=====================

Pan-cancer clinical, gene expression and mutation data {#s4-1}
------------------------------------------------------

The 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 ([@bib47]). 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 ([@bib11]). 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 ([@bib9]). Specifically, GSE39612 ([@bib13]), GSE22396 ([@bib29]), GSE36150 ([@bib24]), GSE50451 ([@bib8]), GSE99316 ([@bib33]) were identified and downloaded.

Implementation of GSVA {#s4-2}
----------------------

APM gene expression status and infiltration levels for immune cell types were quantified using the GSVA method implemented in the R package GSVA ([@bib12]). 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 ([@bib35]). Gene lists for APM score and quantification of immune cell type are provided in [Figure 1---source data 1](#fig1sdata1){ref-type="supplementary-material"} and [Figure 2---source data 1](#fig2sdata1){ref-type="supplementary-material"}.

Calculation of immune infiltration score {#s4-3}
----------------------------------------

The 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 ([@bib27]) and IIS have been described previously ([@bib35]). TIMER ([@bib20]) 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](#fig2s1){ref-type="fig"}).

APM score normalization for TIGS calculation {#s4-4}
--------------------------------------------

Original 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 is$$APS_{normalized} = \ \frac{APS - APS_{pancan\_ min}}{APS_{pancan\_ max} - \ APS_{pancan\_ min}}$$where $APS_{pancan\_ min}$ is the minimal APM score among TCGA pan-cancer samples; and $APS_{pancan\_ max}$ 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.

Normalization of TMB data for TIGS calculation {#s4-5}
----------------------------------------------

TMB was defined as the number of non-synonymous alterations per megabase (Mb) of genome examined. As reported previously ([@bib5]), 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).

TIGS calculation {#s4-6}
----------------

We calculated TIGS as following:$$TIGS = \ APS_{normalized}\  \times \ log\left( {TMB} \right)$$

The 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:$$TIGS = \ APS_{normalized}\  \times \ ln\left( {TMB\  + 1} \right)$$or$$TIGS = \ APS_{normalized}\  \times \ ln\left( {\frac{\ whole\ exome\ mutation\ number}{38} + 1} \right)$$

Immunotherapy clinical studies search strategy {#s4-7}
----------------------------------------------

The dataset search strategy for assessment of cancer immunotherapy ORR) assessment has been described previously ([@bib48]). 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](#fig4sdata1){ref-type="supplementary-material"}. The TMB information for major solid tumor types or subtypes has been described previously ([@bib5]). 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](#fig4sdata2){ref-type="supplementary-material"}). Linear regression models were constructed to correlate ORR with APS, TMB and TIGS for each of the cancer types or subtypes.

Collection and analysis of immunotherapy genomics datasets {#s4-8}
----------------------------------------------------------

To 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 [@bib41] dataset was downloaded from the supplementary files of reference ([@bib41]). 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 ([@bib41]). The [@bib15] dataset was downloaded from the supplementary files of reference ([@bib15]). 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 ([@bib15]). The [@bib38] dataset ([@bib38]) 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 ([@bib38]). 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.

Performance comparison on predicting immunotherapy response {#s4-9}
-----------------------------------------------------------

The 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 [@bib17]. The average expression values among all members defined by the original publications were used to quantify each biomarker. The interferon gamma gene expression signature ([@bib1]) (IFNG) used genes *IFNG*, *STAT1*, *IDO1*, *CXCL10*, *CXCL9*, and *HLA-DRA*. The calculation of IFNG.GS and ISG.RS scores were previously described in [@bib3]. 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.

Statistical analysis {#s4-10}
--------------------

Univariate 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 ([@bib42]). 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 ([@bib32]). 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 ([@bib30]). 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 ([@bib49]) with custom gene sets downloaded from Molecular Signature Database v6.2 ([@bib21]; [@bib39]). 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).

Data availability {#s4-11}
-----------------

All of the code and data used to generate the figures are freely available at <https://github.com/XSLiuLab/tumor-immunogenicity-score> ([@bib44]; 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](#fig1){ref-type="fig"}, [2](#fig2){ref-type="fig"}, [4](#fig4){ref-type="fig"} and [5](#fig5){ref-type="fig"}.

Funding Information
===================

This paper was supported by the following grant:

-   http://dx.doi.org/10.13039/501100001809National Natural Science Foundation of China 31771373 to Xue-Song Liu.

We thank the authors and participating patients of the immunotherapy publications for providing the data used for this analysis. Our gratitude is also extended to the TCGA project for making cancer genomics data available for analysis. We thank Raymond Shuter for editing the text. Thanks also to the ShanghaiTech University High Performance Computing Public Service Platform for providing computing services. Thanks also to other members of Liu lab for helpful discussions.

Additional information {#s5}
======================

No competing interests declared.

Resources, Data curation, Software, Formal analysis, Validation, Investigation, Visualization, Methodology.

Investigation, Methodology.

Investigation, Methodology.

Investigation, Methodology.

Conceptualization, Resources, Formal analysis, Supervision, Funding acquisition, Validation, Investigation, Methodology, Writing---original draft, Project administration, Writing---review and editing.

Additional files {#s6}
================

10.7554/eLife.49020.022

Data availability {#s7}
=================

All 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.

The following previously published datasets were used:

HarmsPBichakjianC2013Distinct gene expression profiles of viral- and non-viral associated Merkel cell carcinoma revealed by transcriptome analysisNCBI Gene Expression OmnibusGSE39612

PaulsonKGIyerJGSchelterJClearyMAHardwickJNghiemP2011Gene expression analysis of Merkel Cell CarcinomaNCBI Gene Expression OmnibusGSE22396

MastersonLThibodeauBJFortierLEGeddesTJPruetzBLKeidanRWilsonGD2014Gene expression changes associated with prognosis of Merkel cell carcinomaNCBI Gene Expression OmnibusGSE36150

BrownellIDailyK2015Microarray analysis of Merkel cell carcinoma (MCC) tumors, small cell lung cancer (SCLC) tumors, and MCC cell linesNCBI Gene Expression OmnibusGSE50451

SatoTKanedaATsujiSIsagawaTYamamotoSFujitaTYamanakaRTanakaYNukiwaTMarquezVEIshikawaYIchinoseMAburataniH2013Gene repression and ChIP-seq in Human Small Cell Lung CancerNCBI Gene Expression OmnibusGSE99316

10.7554/eLife.49020.035

Decision letter

Freeman

Gordon

Reviewing Editor

Dana-Farber Cancer Institute

United States

Ji

Hongbin

Reviewer

Institute of Biochemistry and Cell Biology, Shanghai Institutes for Biological Sciences, Chinese Academy of Sciences

China

In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.

**Acceptance summary:**

Identification of biomarkers that can predict which patients will benefit from immune checkpoint inhibition therapy is clinically important. Wang et al., describe a new computational method to identify responders to immune checkpoint inhibitors by calculating a tumor immunogenicity score (TIGS). TIGS combines tumor mutational burden (TMB) with a gene set score of 18 genes associated with MHC class I Antigen Presentation Machinery (APM) score. They describe the APM score across cancer types in TCGA and correlate APM with other gene expression pathways and immune cell infiltration across cancers. In both a pan-cancer analysis of ICI objective response rates and an ICI clinical response prediction for individual patients, they show that TIGS predicts response to ICI better than TMB alone, PD-L1, immune infiltrate, or Interferon-gene signatures, and somewhat better than the TIDE method based on T cell dysfunction and exclusion gene expression signature. TIGS is a tumor inherent biomarker and may be valuable in predicting response to immunotherapy as well as guiding ways to enhance the immunogenicity of tumors.

**Decision letter after peer review:**

Thank you for submitting your article \"Antigen presentation and tumor immunogenicity in cancer immunotherapy response prediction\" for consideration by *eLife*. Your article has been reviewed by three peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Tadatsugu Taniguchi as the Senior Editor. The following individual involved in review of your submission has agreed to reveal their identity: Hongbin Ji (Reviewer \#1).

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

Summary:

Wang et al. have developed a method named tumor immunogenicity score (TIGS) that combines tumor mutational burden (TMB) and antigen processing and presenting machinery gene expression signature to measure tumor immunogenicity. They found that TIGS could outperform TMB and other known ICI response prediction biomarkers in both pan-cancer ICI objective response rates correlation and ICI clinical response prediction.

Essential revisions:

1\) The results of Şenbabaoğlu should be cited more fully and their 18 gene APS compared to the 7 gene Şenbabaoğlu APM geneset. The antigen presenting set is all MHC class I. What is the reasoning behind not including MHC class II genes like DRB1, DRB2, CIITA?

2\) TIGS should be compared to TIDE + TMB. Inflamed gene expression signatures have not been shown to be a highly predictive biomarker so I think that comparisons to interferon signatures will show superiority but not guide any advance in our thinking.

3\) Multiple reviewers were perplexed by why high TIGs is associated with good outcome in some tumor types but poorer outcomes in other tumor types (Figure 3B). Expand the discussion here.

4\) Improve the Cox regression analysis as suggested and make the definition of high or low marker expression consistent throughout.

5\) The suggestion to apply the new analyses to many cancer types for ICI response prediction is limited by the available datasets that have all of the required information. The 3 datasets analyzed are what is available and these have been done.

6\) Indicate the source or reference for their linear correlation formula -- \"objective response rate = 21.4 ×TIGS -- 2.7, \".

7\) If a patient has a beta2m mutation, is this captured by the APM signature and TIGS method?

Reviewer \#1:

Identification of the biomarkers that predict which patients may benefit from immune checkpoint inhibition therapy is clinically important. In this study, Wang et al. developed a method named tumor immunogenicity score (TIGS) that combined tumor mutational burden (TMB) and antigen processing and presenting machinery gene expression signature to measure tumor immunogenicity. They found that TIGS could outperform TMB and other known ICI response prediction biomarkers in both pan-cancer ICI objective response rates correlation and ICI clinical response prediction. Thus, they proposed that TIGS is a potential tumor inherent biomarker for ICI response prediction. Overall, the study is novel and interesting. I list some of my concerns below.

1\) The definition of high or low marker expression should be consistent throughout. For example, the authors defined patients with APS of first quartile as \"APS-High\", and those at the fourth quartile as \"APS-Low\". In contrast, they defined patients with TIGS above the median as \"TIGS-high\", and the remaining as \"TIGS-low\". Similar issues exist for the definition of TMB. Please make the correction and explain the rationale behind.

2\) In Figure 3B, Cox regression analysis show that high TIGS is significantly associated with poor survival of patients in several types of malignancies such as adrenocortical carcinoma (HR=5.23, p=0.00105), Kidney Chromophobe (HR=89.9, p=0.01408), Thymoma (HR=8.22, p=0.00198)...etc. How to explain this phenomena given that TIGS reflects tumor immunogenicity and high TIGS predicts favorable prognosis in patients following immunotherapy.

3\) In this study, the authors demonstrate an improved predictive power of TIGS in ICI clinical response when compared to TMB and other gene expression profiling based biomarkers, such as TIDE, IIS, IFNΓ. The authors should discuss the potential mechanisms underlying the superior performance of TIGS for immunotherapy clinical response prediction. The authors state \"Furthermore, our linear correlation formula -- objective response rate = 21.4 ×TIGS -- 2.7, -- can be used to....\". Please also indicate the source or reference for this formula.

Reviewer \#2:

This well-written article describes a score, the Tumor Immunogenicity Score (TIGS) which combines tumor mutation burden (TMB) and a gene set score of 18 genes associated with Antigen Presentation Machinery (APM) score (APS). The 18 genes were described by Leone et al., 2013) and are PSMB5, PSMB6, PSMB7, PSMB8, PSMB9, PSMB10, TAP1, TAP2, ERAP1, ERAP2, CANX, CALR, PDIA3, TAPBP, B2M, HLA-A, HLA-B, and HLA-C. The author show that the TIGS (product of the natural log of TMB and the normalized APS) marginally outperforms TMB or APS in A) linear associated with Objective Response Rate (ORR), which is a rank order of tissue response to PD1 and PDL1 therapy and B) prediction of response in 2 tissue (Melanoma, Urothelial) in 3 studies (Van Allen, 2015, Hugo et al., 2016 and Snyder et al., 2017).

This article represents an extension to Şenbabaoğlu et al., 2016 which showed that a 7-gene APM signature of (HLA-A/B/C, B2M, TAP1, TAP2, and TAPBP), that is contained in the authors 18-gene signature was associated of immune infiltration, TMB, and ORR to immune therapy in kidney. Many of the figures presented by the author duplicate or extend that article.

Some of the figures of APM and APS are redundant with Şenbabaoğlu et al., and I would suggest that these should be moved to the supplement. Instead could the authors should provide a more in-depth discussion of the TIGS scores, IFN etc. Figure 1 should be moved to the supplement.

Are the APS results different to the 7-gene signature described by Şenbabaoğlu et al., There should be a comparison of the performance of the 7-gene Şenbabaoğlu and the 18-gene signature here.

Which cells (CIBERSORT scores) or subsets of the IIS scores are most associated with the APS scores.

The authors should compare their TIGS/APS to recent scores for immune presentation PHBR scores for MHC I and II (Marty et al., 2017 and 2018). https://www.cell.com/cell/pdf/S0092-8674(18)31109-7.pdf https://www.cell.com/cell/pdf/S0092-8674(17)31144-3.pdf

Does APS correlate with PHBR I alone or does it capture both PHBR scores for MHC class I and II? Does TIGS outperform a product of PHBR and TMB?

The authors recently stated that TMB is associated with gender (Wang et al., 2019b). What is the association between TIGS, APS and gender.

In the Introduction, some of the statements are over-generalized and do not reflect the complexity in defining good predictors of immunotherapy response. Whilst a correlation exists, TMB does not always predict response, neither does TIL. Some cancers (e.g. renal) with high immune infiltrate have poor response. Please edit the Introduction to reduce broad over simplifications or generalizations.

Reviewer \#3:

Wang et al., describe a new computational method to identify responders for immune checkpoint inhibitors (ICI) using gene expression data from the TCGA database. The authors calculate a tumor immunogenicity score (TIGS) by combining tumor mutational burden (TMB) with antigen processing and presentation machinery (APM) gene signatures. They describe APM signatures across cancer types in TCGA and correlate APM with other gene expression pathways and immune cell infiltration across cancers. They next evaluate the ability of TIGS to predict response to ICI and show improved predictions using this method, compared to TMB alone, the TIDE method by Liu and colleagues, and several biomarkers (PDL1, CD8, etc.). This study is timely given the broad interest in predicting clinical responses to ICI therapy, and the concept of combining antigen presentation gene expression with TMB is also novel and interesting. However, the study is lacking in benchmarking data that support the use of the APM gene signature, and in comprehensive comparisons to prior gene signatures that also synergize with TMB in predicting immune response to ICI. It is also unclear whether the authors have evaluated the performance of prior gene signatures combined with TMB, compared to the TIGS method. Without these comparisons, it is difficult for the reader to truly evaluate the value added by this method, and I suspect will result in lower adoption of the method by the community.

1\) A major premise of using antigen presentation gene scores is the presence of mutations in a subset of these genes in non-responder patients (i.e. Zaretsky et al., 2016). Therefore, the authors should determine whether their APM analysis is able to actually capture these defects in tumor samples. In other words, if a patient has a b2m mutation, is this captured by the APM signature and TIGS method? Are these instances the major driver of the value of APM analysis, or are changes in expression levels (without mutation or LOH) also predictive of response?

2\) The authors have compared the performance of TIGS to several other prediction tools, however this analysis is not fully described, and I have several questions:

\- The comparisons to TIDE are interesting. As I understand it, TIDE only takes into account gene expression, and not TMB. In contrast, TIGS takes into account gene expression (APM) and TMB. The authors should show the data for APM signature alone in several of the figures, for example in Figure 5A-C.

\- The authors should also clarify in the main text whether TIDE incorporates TMB into their calculations, and if not, the authors should compare the performance of TIDE +TMB to TIGS.

\- Similarly, for PDL1, IFNΓ, and CD8 scores, were these also combined with TMB? Or were they used in isolation to predict response (Figure 5)? The question is: what is the real value added -- is it the APM score, or combining gene expression with TMB?

\- If TIGS remains a better predictor of response compared to TIDE + TMB, the authors should describe in a main figure the performance comparison in all TCGA cancer types, rather than showing the comparison in 3 (2 that currently perform similarly, and 1 where TIGS outperforms). This information, in a main figure, is critical for the reader to understand the value of this new method across many cancers.

\- Since APM genes are turned on by the IFN pathway (as the authors discuss), I would like to see more comprehensive comparisons to IFN pathway signature predictions, beyond only the 6 IFNΓ gene signature score taken from the Ayers et al. manuscript. In particular, I would like to see comparisons to the ISG.RS and IFNΓ.GS signatures described in Benci et al., 2019.

3\) The authors correlate APM score with immune infiltration (using IIS and TIGER). Given the prior concerns regarding data normalization in the TIGER method (Newman et al., 2017), I would suggest adding an additional comparison using CIBERSORT. The correlation of APM and immune cell infiltration is independently interesting (without the prediction of response rates), and I think it would be useful to dig into this a bit more -- i.e. which cell types correlate most with high APM scores?

10.7554/eLife.49020.036

Author response

> Essential revisions:
>
> 1\) The results of Şenbabaoğlu should be cited more fully and their 18 gene APS compared to the 7 gene Şenbabaoğlu APM geneset. The antigen presenting set is all MHC class I. What is the reasoning behind not including MHC class II genes like DRB1, DRB2, CIITA?

As pointed out by the second reviewer, Şenbabaoğlu et al., 2016 performed antigen presentation gene expression signature analysis, our method for analyzing antigen presentation gene expression signature is similar to Şenbabaoğlu et al., 2016 but with different gene list. And this new information about the citation of Şenbabaoğlu et al., 2016 have been included in the first section of Results in the revised manuscript (subsection "APM score definition and pan-cancer analysis", last paragraph).

In addition, as suggested, we compared the performance of 7 genes in Şenbabaoğlu et al., 2016 with our 18 genes, results suggested that these two methods show strong association in TCGA pan-cancer level (new Figure 1---figure supplement 1) and similar performance in ICI response prediction ([Author response image 1](#respfig1){ref-type="fig"}).

![](elife-49020-resp-fig1){#respfig1}

MHC I are found on the cell surface of all nucleated cells, and function in displaying peptide fragments of proteins from within the cell to cytotoxic CD8^+^ T cells. MHC II are normally found only on professional antigen-presenting cells such as dendritic cells. The antigens presented by MHC II are derived from extracellular proteins. The anti-cancer immune response against mutated peptides (neoantigen) is primarily attributed to MHC-I-restricted cytotoxic CD8^+^ T cell responses. MHC-II-restricted CD4^+^ T cells also drive anti-tumor responses, however their contribution to neoantigen presentation is not clear.

This manuscript primarily focused on the cytosolic or endogenous neoantigen presentation pathway mediated by MHC I, this does not mean that the potential neoantigen presentation by MHC II is not important. We agree with the reviewers that MHC II presentation could also contribute to the immunogenicity of cancer cells. And the related new discussion has also been included in the revised Discussion part (Discussion, fifth paragraph).

> 2\) TIGS should be compared to TIDE + TMB. Inflamed gene expression signatures have not been shown to be a highly predictive biomarker so I think that comparisons to interferon signatures will show superiority but not guide any advance in our thinking.

TIDE reflect gene expression signatures of T cell dysfunction and exclusion, it may do not have a similar rationale to be combined with TMB as APS in this study. As suggested, we now include the combination of TIDE and TMB in ICI clinical response prediction ([Author response image 1](#respfig1){ref-type="fig"}). In the Snyder et al., 2017 dataset, combination of TIDE and TMB still show poor ICI response prediction.

> 3\) Multiple reviewers were perplexed by why high TIGs is associated with good outcome in some tumor types but poorer outcomes in other tumor types (Figure 3B). Expand the discussion here.

Most TCGA cancer patients have not been treated with immunotherapy, and the prognosis of these cancer patients are influenced by many factors. Different cancer types could have different prognosis in regards to high APS, TMB or TIGS.

Meta-analysis with all TCGA cancer types suggests that APS is not associated with cancer patients' prognosis (new Figure 1B), patients with high TMB tends to have poor prognosis (new Figure 2---figure supplement 3), and this observation is similar to previous studies (Owada-Ozaki et al. Prognostic Impact of Tumor Mutation Burden in Patients With Completely Resected Non--Small Cell Lung Cancer: Brief Report. J Thorac Oncol. 2018 Aug;13(8):1217-1221; McNamara et al. Prognostic and predictive impact of high tumor mutation burden (TMB) in solid tumors: A systematic review and meta-analysis. Annals of Oncology 30 (Supplement 5): v25--v54, 2019).

TMB reflect tumor antigenicity, also predict improved survival after immunotherapy. However in cancer patients not treated with immunotherapy, high TMB tends to be associated with poor prognosis, probably because tumor accumulate mutation during progression due to genome instability, and consequently high TMB is usually associated with late stage of cancer.

High TIGS also tends to associated with poor prognosis. This new data and discussion have now been included in the revised manuscript (new Figure 3B). The poor prognosis associated with high TIGS in cancer patients not treated with immunotherapy may be due to similar mechanism as high TMB. This new data analysis and discussion have now been included in the revised manuscript (subsection "Tumor immunogenicity score (TIGS) definition and pan-cancer profiling", last paragraph).

> 4\) Improve the Cox regression analysis as suggested and make the definition of high or low marker expression consistent throughout.

Cox regression analysis has now been improved as suggested (see revised Figure 1B, Figure 2---figure supplement 3 and Figure 3B). The definition of high and low markers have now been consistent, and Figure 2 and Figure 2---figure supplement 1 has been edited based on new definition.

> 5\) The suggestion to apply the new analyses to many cancer types for ICI response prediction is limited by the available datasets that have all of the required information. The 3 datasets analyzed are what is available and these have been done.

Thanks for this point, we replied the third reviewer as suggested.

> 6\) Indicate the source or reference for their linear correlation formula -- \"objective response rate = 21.4 ×TIGS -- 2.7,\".

This is based on the linear regression model of Figure 4, and this information has now been included in the revised manuscript (Discussion, third paragraph).

> 7\) If a patient has a beta2m mutation, is this captured by the APM signature and TIGS method?

Our original method for calculating APM score was based on mRNA expression, it does not directly capture mutation status of APM genes. If a patient has a loss of function beta2m mutation, theoretically, this patient will lose the ability to present antigens through MHC I to immune system, and in the revised APS quantification method, the APS will be reset to zero. Loss of function mutation in beta2m is very rare, and we could not find it in both TCGA and three ICI datasets. Therefore changes in expression level appear to be major driver for APS differences.

> Reviewer \#1:
>
> \[...\] 1) The definition of high or low marker expression should be consistent throughout. For example, the authors defined patients with APS of first quartile as \"APS-High\", and those at the fourth quartile as \"APS-Low\". In contrast, they defined patients with TIGS above the median as \"TIGS-high\", and the remaining as \"TIGS-low\". Similar issues exist for the definition of TMB. Please make the correction and explain the rationale behind.

We thank the reviewer for this point. As suggested, we now use the median value as cutoff throughout this study, patients with APS or TMB or TIGS values above the median value were defined as APS-high or TMB-high or TIGS-high respectively. The Figure 2 and Figure 2---figure supplement 1 have thus been corrected based on this new cutoff. Generally the new data do not show apparent difference compared with original one.

> 2\) In Figure 3B, Cox regression analysis show that high TIGS is significantly associated with poor survival of patients in several types of malignancies such as adrenocortical carcinoma (HR=5.23, p=0.00105), Kidney Chromophobe (HR=89.9, p=0.01408), Thymoma (HR=8.22, p=0.00198)...etc. How to explain this phenomena given that TIGS reflects tumor immunogenicity and high TIGS predicts favorable prognosis in patients following immunotherapy.

We thank the reviewer for this point. Most TCGA cancer patients have not been treated with immunotherapy, and the prognosis of these cancer patients are influenced by many factors. Different cancer types could have different prognosis in regards to high APS, TMB or TIGS.

Meta-analysis with all TCGA cancer types suggests that APS is not associated with cancer patients' prognosis (new Figure 1B), patients with high TMB tends to have poor prognosis (new Figure 2---figure supplement 3), and this observation is similar to previous studies (Owada-Ozaki et al. Prognostic Impact of Tumor Mutation Burden in Patients With Completely Resected Non--Small Cell Lung Cancer: Brief Report. J Thorac Oncol. 2018 Aug;13(8):1217-1221; McNamara et al. Prognostic and predictive impact of high tumor mutation burden (TMB) in solid tumors: A systematic review and meta-analysis. Annals of Oncology 30 (Supplement 5): v25--v54, 2019).

TMB reflect tumor antigenicity, also predict improved survival after immunotherapy. However in cancer patients not treated with immunotherapy, high TMB tends to be associated with poor prognosis, probably because tumor accumulate mutations during progression due to genome instability, and consequently high TMB is usually associated with late stage of cancer.

High TIGS also tends to associated with poor prognosis. This new data and discussion have now been included in the revised manuscript (new Figure 4B). The poor prognosis associated with high TIGS in cancer patients not treated with immunotherapy may be due to similar mechanism as high TMB. This new data analysis and discussion have now been included in the revised manuscript (subsection "Tumor immunogenicity score (TIGS) definition and pan-cancer profiling", last paragraph).

> 3\) In this study, the authors demonstrate an improved predictive power of TIGS in ICI clinical response when compared to TMB and other gene expression profiling based biomarkers, such as TIDE, IIS, IFNΓ. The authors should discuss the potential mechanisms underlying the superior performance of TIGS for immunotherapy clinical response prediction. The authors state \"Furthermore, our linear correlation formula -- objective response rate = 21.4 ×TIGS -- 2.7, -- can be used to....\". Please also indicate the source or reference for this formula.

TIGS captured the two key aspects of tumor immunogenicity, tumor antigenicity and antigen presentation, and theoretically this tumor immunogenicity score should have improved prediction power compared to biomarkers that only reflect tumor antigenicity (such as TMB), or tumor immune environment (such as TIDE IFNΓ). As suggested, this new discussion has now been included (Discussion, first paragraph). In the Results section, we also describe the rationale of TIGS (subsection "Tumor immunogenicity score (TIGS) definition and pan-cancer profiling"). The linear correlation formula was based on data in Figure 4C, and this information has been included in the Discussion part (third paragraph).

> Reviewer \#2:
>
> This well-written article describes a score, the Tumor Immunogenicity Score (TIGS) which combines tumor mutation burden (TMB) and a gene set score of 18 genes associated with Antigen Presentation Machinery (APM) score (APS). The 18 genes were described by Leone et al., 2013) and are PSMB5, PSMB6, PSMB7, PSMB8, PSMB9, PSMB10, TAP1, TAP2, ERAP1, ERAP2, CANX, CALR, PDIA3, TAPBP, B2M, HLA-A, HLA-B, and HLA-C. The author show that the TIGS (product of the natural log of TMB and the normalized APS) marginally outperforms TMB or APS in A) linear associated with Objective Response Rate (ORR), which is a rank order of tissue response to PD1 and PDL1 therapy and B) prediction of response in 2 tissue (Melanoma, Urothelial) in 3 studies (Van Allen, 2015, Hugo et al., 2016 and Snyder et al., 2017).
>
> This article represents an extension to Şenbabaoğlu et al., 2016 which showed that a 7-gene APM signature of (HLA-A/B/C, B2M, TAP1, TAP2, and TAPBP), that is contained in the authors 18-gene signature was associated of immune infiltration, TMB, and ORR to immune therapy in kidney. Many of the figures presented by the author duplicate or extend that article.

We agree with the reviewer that Şenbabaoğlu et al., 2016 performed a comprehensive study about cancer immune cell infiltration. Şenbabaoğlu et al., 2016 also reported that elevated 7-APM gene expression signature could predict clinical response to Nivolumab (anti-PD-1) in a ccRCC clinical trial, however, this clinical trial only has 6 patients, and their conclusion need to be validated with more patients. Our study focused the immunogenicity of tumor cells, and propose to combine antigen presentation and tumor mutation burden together in quantifying tumor immunogenicity, and thus the major ideas and conclusions of our study are different from Şenbabaoğlu et al., 2016

> Some of the figures of APM and APS are redundant with Şenbabaoğlu et al., 2016 and I would suggest that these should be moved to the supplement. Instead could the authors should provide a more in-depth discussion of the TIGS scores, IFN etc. Figure 1 should be moved to the supplement.

We prefer to keep Figure 1 as main figure, Şenbabaoğlu et al., 2016 did show the comparison of APM between normal vs. cancer in 15 TCGA cancer types in their Supplementary Figure 11. The methods for APM calculation have some similarity, however, the focus and data of our Figure 1 and their Supplementary Figure 11 are different. Their study focused on the comparison of APM between cancer vs. normal, however our study focused on the pan-cancer distribution of APS, and this is critical for further compare the pan-cancer ORR in Figure 4. We include the citation of Şenbabaoğlu et al.'s Supplementary Figure 11 data, which show the comparison of APM between cancer vs. normal, in the Results section (subsection "APM score definition and pan-cancer analysis, last paragraph).

> Are the APS results different to the 7-gene signature described by Şenbabaoğlu et al., 2016 There should be a comparison of the performance of the 7-gene Şenbabaoğlu and the 18-gene signature here.

The APS generated with our 18 genes and with 7-gene signature described by Şenbabaoğlu et al., 2016 are highly correlated (Figure 1---figure supplement 1), they also show similar prediction power in immunotherapy response prediction ([Author response image 1](#respfig1){ref-type="fig"}).

> Which cells (CIBERSORT scores) or subsets of the IIS scores are most associated with the APS scores.

We thank the reviewer for this point. As suggested we now included the correlation analysis between APS and cell type status calculated with both CIBERSORT and IIS (new Figure 2---figure supplement 4).

> The authors should compare their TIGS/APS to recent scores for immune presentation PHBR scores for MHC I and II (Marty et al., 2017 and 2018). https://www.cell.com/cell/pdf/S0092-8674(18)31109-7.pdf https://www.cell.com/cell/pdf/S0092-8674(17)31144-3.pdf
>
> Does APS correlate with PHBR I alone or does it capture both PHBR scores for MHC class I and II? Does TIGS outperform a product of PHBR and TMB?

As suggested, we included the comparison between APS and PHBR in the revised manuscript (new Figure 1---figure supplement 1).

Our APS captures the expression level information of MHC genes in patient level. Patient Harmonic Best Rank (PHBR) score represents antigen presentation ability for mutations. Both PHBR I and II scores are determined by the MHC genotypes of patients. We obtained patient-by-mutation PHBR-I/II score matrix for TCGA patients from the author of PHBR-I/II paper and summarized the median for each patient, followed by calculation of the correlation between APS and PHBR sores in pan-cancer level. The result shows that APS do not correlate with both PHBR I and PHBR II scores (new Figure 1---figure supplement 1).

PHBR could not capture the antigen presentation differences that exist in different cancer types, since different cancer types should have exactly the same PHBR status because PHBR are MHC genotypes based, thus PHBR cannot be used to explain the ORR differences in cancer types of different tissue origin. On the contrary, our APS and TIGS are gene expression based, and can be used in ORR prediction in different cancer types (Figure 4).

For ICI clinical response prediction in individual cancer patients, raw sequencing information is not available, and this is required for both HLA typing and mutation calling, thus we reviewed the supplementary files of three immunotherapy datasets and found that only two of them (Van Allen, 2015 and Snyder, 2017) have both MHC-I and mutation information available. We then calculated the PHBR I score for each residue and summarized the median for each patient. The result score was used to predict the immunotherapy response. The analysis process is recorded in https://github.com/XSLiuLab/pypresent/tree/master/icb_analysis. The result shows that TIGS outperform PHBR I alone or PHBR I combined with TMB in immunotherapy clinical response prediction ([Author response image 1](#respfig1){ref-type="fig"}).

> The authors recently that TMB is associated with gender (Wang et al., 2019b). What is the association between TIGS, APS and gender.

In our IJC 2019 paper, we did observe a significant gender difference in TMB's prediction power in lung cancer (349 patient, 171 men 178 women). Due to lack of sufficient number of patients, we could not draw significant gender difference in other cancer types. The three ICI datasets with both genomic DNA mutation and RNA expression data available are melanoma and urothelial cancer. In these datasets, we do not have sufficient number of patients to investigate the gender difference in biomarker's performance. Related discussion has been included in the revised manuscript (Discussion, fifth paragraph).

> In the Introduction, some of the statements are over-generalized and do not reflect the complexity in defining good predictors of immunotherapy response. Whilst a correlation exists, TMB does not always predict response, neither does TIL. Some cancers (e.g. renal) with high immune infiltrate have poor response. Please edit the Introduction to reduce broad over simplifications or generalizations.

We thank the reviewer for this point. We edit the Introduction part as suggested (first paragraph).

> Reviewer \#3:
>
> Wang et al., describe a new computational method to identify responders for immune checkpoint inhibitors (ICI) using gene expression data from the TCGA database. The authors calculate a tumor immunogenicity score (TIGS) by combining tumor mutational burden (TMB) with antigen processing and presentation machinery (APM) gene signatures. They describe APM signatures across cancer types in TCGA and correlate APM with other gene expression pathways and immune cell infiltration across cancers. They next evaluate the ability of TIGS to predict response to ICI and show improved predictions using this method, compared to TMB alone, the TIDE method by Liu and colleagues, and several biomarkers (PDL1, CD8, etc.). This study is timely given the broad interest in predicting clinical responses to ICI therapy, and the concept of combining antigen presentation gene expression with TMB is also novel and interesting. However, the study is lacking in benchmarking data that support the use of the APM gene signature, and in comprehensive comparisons to prior gene signatures that also synergize with TMB in predicting immune response to ICI. It is also unclear whether the authors have evaluated the performance of prior gene signatures combined with TMB, compared to the TIGS method. Without these comparisons, it is difficult for the reader to truly evaluate the value added by this method, and I suspect will result in lower adoption of the method by the community.
>
> 1\) A major premise of using antigen presentation gene scores is the presence of mutations in a subset of these genes in non-responder patients (i.e. Zaretsky et al., 2016). Therefore, the authors should determine whether their APM analysis is able to actually capture these defects in tumor samples. In other words, if a patient has a b2m mutation, is this captured by the APM signature and TIGS method? Are these instances the major driver of the value of APM analysis, or are changes in expression levels (without mutation or LOH) also predictive of response?

We thank the reviewer for this point. Our original method for calculating APM score was based on mRNA expression, it does not directly capture mutation status of APM genes. If a patient has a loss of function b2m mutation, theoretically, this patient will lose the ability to present antigens to immune system, and in the revised APM quantification method, the APS will be reset to zero. Loss of function mutation in b2m is very rare, and we could not find it in both TCGA and three ICI datasets. Therefore changes in expression level appear to be major driver for APM differences.

> 2\) The authors have compared the performance of TIGS to several other prediction tools, however this analysis is not fully described, and I have several questions:
>
> \- The comparisons to TIDE are interesting. As I understand it, TIDE only takes into account gene expression, and not TMB. In contrast, TIGS takes into account gene expression (APM) and TMB. The authors should show the data for APM signature alone in several of the figures, for example in Figure 5A-C.

We actually showed the performance of APS in ICI response prediction in the original Figure 5D/E/F.

> \- The authors should also clarify in the main text whether TIDE incorporates TMB into their calculations, and if not, the authors should compare the performance of TIDE +TMB to TIGS.

TIDE reflect gene expression signatures of T cell dysfunction and exclusion, it may be do not have a similar rationale to be combined with TMB as APS in this study. As suggested we combined TIDE with TMB, and this combination still shows poor predictive power in Snyder et al., 2017 dataset ([Author response image 1](#respfig1){ref-type="fig"}).

> \- Similarly, for PDL1, IFNΓ, and CD8 scores, were these also combined with TMB? Or were they used in isolation to predict response (Figure 5)? The question is: what is the real value added -- is it the APM score, or combining gene expression with TMB?

The combination between APS and TMB has been driven by the rationale that tumor immunogenicity can be divided into two independent steps: neoantigen generation through mutation and antigen presentation through MHC, and this information has been described in the Introduction and Results section. In the original manuscript PDL1, IFNΓ and CD8 score has not been combined with TMB, since we could not find proper rationale for the suggested combination between PDL1, IFNΓ or CD8 with TMB. The key information in Figure 5 is that TIGS measured with combined APS and TMB can outperform known biomarkers in ICI clinical response prediction.

> \- If TIGS remains a better predictor of response compared to TIDE + TMB, the authors should describe in a main figure the performance comparison in all TCGA cancer types, rather than showing the comparison in 3 (2 that currently perform similarly, and 1 where TIGS outperforms). This information, in a main figure, is critical for the reader to understand the value of this new method across many cancers.

Currently we only have three datasets available for analysis. Since ICI datasets that include patients' genomic data, gene mRNA expression data and clinical response data are very limited, currently only these three datasets are available.

> \- Since APM genes are turned on by the IFN pathway (as the authors discuss), I would like to see more comprehensive comparisons to IFN pathway signature predictions, beyond only the 6 IFNΓ gene signature score taken from the Ayers et al. manuscript. In particular, I would like to see comparisons to the ISG.RS and IFNΓ.GS signatures described in Benci et al., 2019.

As suggested, we now included ISG.RS and IFNΓ.GS in ICI clinical response prediction comparison (new Figure 5D/E/F).

> 3\) The authors correlate APM score with immune infiltration (using IIS and TIGER). Given the prior concerns regarding data normalization in the TIGER method (Newman et al., 2017), I would suggest adding an additional comparison using CIBERSORT. The correlation of APM and immune cell infiltration is independently interesting (without the prediction of response rates), and I think it would be useful to dig into this a bit more -- i.e. which cell types correlate most with high APM scores?

We thank the reviewer for this point. As suggested, we now include CIBERSORT analysis (new Figure 2---figure supplement 4). The cell types correlated with high APS have now been investigated through both CIBERSORT and IIS analysis (new Figure 2---figure supplement 4).
