Integrated Genomic Characterization of the Human Immunome in Cancer Free

Immunotherapy has been considered as a promising strategy for treatment of various types of cancer. Although these immunotherapy methods have made remarkable success, only a minority of patients are observed to respond to the treatment with cytotoxic T lymphocyte antigen-4 (CTLA-4) or programmed death receptor-1 (PD-1) blockade (1). There is a clinical need to identify predictive biomarkers and to understand the potential mechanisms of immunotherapy resistance.

With the development of high-throughput sequencing technology, recent evidence has pointed to a number of predictive biomarkers, such as tumor mutational burden (TMB; ref. 2), neoantigen burden (3), PD-L1 or PD-L2 mRNA expression (4), epigenetic markers, and immune cell infiltration profiles (5). With these predictive biomarkers in hand, machine learning–based methods are being built to predict the immunotherapy response. However, it is challenging to integrate different biomarkers. Moreover, increasing studies are beginning to reveal the underlying mechanisms of immunotherapy resistance. The expression of immune-related genes (6), insufficient immune cell infiltration (7), oncogenic pathway activity (8) as well as metabolism dysregulation (9) have all been related to therapy resistance. Nonetheless, additional insights are still needed for a comprehensive understanding of the mechanisms of immunotherapy resistance.

In addition, genes do not function in isolation but rather interact extensively with each other, and therefore, any genetic abnormality is not restricted to the gene product that it encodes (10). The emerging tools of network- or pathway-based medicine offer a valuable platform to explore the predictive biomarkers as well as to understand the dysregulation pattern of cancer-related pathways. Previous studies by The Cancer Genome Atlas (TCGA) have systematically analyzed the alteration landscape of cancer-related signaling pathways (11). However, we still lack the knowledge about the alteration pattern of immune-related pathways across cancer types.

To further refine both the genomic and transcriptomic-derived events associated with immune infiltration and response to immunotherapy, we first charted the detailed perturbation landscape of immune-related pathways. A network-based immunogenomics model with scoring systems was designed to identify modulation in signaling networks associated with immunogenicity. Furthermore, using a network propagation framework integrated with drug-associated gene signatures, we discovered potential immunomodulatory drug candidates. These results can help elucidate the functionally relevant mechanisms of immune-related pathway alterations and might inform potential immunotherapy treatment options.

Human immune-related genes were obtained from the ImmPort project (http://www.immport.org; ref. 12). In total, 1,811 genes in 17 immune-related pathways were obtained. In addition, we also obtained genes that are identified to be essential to immune response from recent literature (13, 14). These genes were defined as essential immunotherapy genes and grouped into the 18th pathway. In total, 2,273 immune-related genes were obtained.

The results in our analysis are based upon datasets generated by TCGA Research Network (http://cancergenome.nih.gov/). We analyzed 33 different TCGA projects; each project represents a specific cancer type, including KIRC, kidney renal clear cell carcinoma; KIRP, kidney renal papillary cell carcinoma; KICH, kidney chromophobe; LGG, brain lower grade glioma; GBM, glioblastoma multiforme; BRCA, breast cancer; LUSC, lung squamous cell carcinoma; LUAD, lung adenocarcinoma; READ, rectum adenocarcinoma; COAD, colon adenocarcinoma; UCS, uterine carcinosarcoma; UCEC, uterine corpus endometrial carcinoma; OV, ovarian serous cystadenocarcinoma; HNSC, head and neck squamous carcinoma; THCA, thyroid carcinoma; PRAD, prostate adenocarcinoma; STAD, stomach adenocarcinoma; SKCM, skin cutaneous melanoma; BLCA, bladder urothelial carcinoma; LIHC, liver hepatocellular carcinoma; CESC, cervical squamous cell carcinoma and endocervical adenocarcinoma; ACC, adrenocortical carcinoma; PCPG, pheochromocytoma and paraganglioma; SARC, sarcoma; LAML, acute myeloid leukemia; PAAD, pancreatic adenocarcinoma; ESCA, esophageal carcinoma; TGCT, testicular germ cell tumors; THYM, thymoma; MESO, mesothelioma; UVM, uveal melanoma; DLBC, lymphoid neoplasm diffuse large B-cell lymphoma; CHOL, cholangiocarcinoma.

Somatic mutations were obtained from the publicly available TCGA MAF file (“MC3”), which covers 10,224 patients (15). This dataset was directly downloaded from Synapse under the number of syn7214402. Six calling methods were applied and numbers of filters were applied. All these mutations were subjected to ANNOVAR (16) and we obtained the conservation and MetaSVM score for each mutation (17).

RNA sequencing (RNA-seq) data were obtained from the TCGA project via the R-package “TCGAbiolinks,” which was specifically developed for integrative analysis with GDC data. We downloaded the Fragments Per Kilobase of transcript per Million mapped reads (FPKM)–based gene expression for 33 types of cancer. In addition, we also obtained the raw read counts for each gene in these samples. The clinical information for patients of 33 cancer types was downloaded from TCGA project via the R-package “TCGAbiolinks,” including the survival status, stages, grades, and survival time.

Moreover, we obtained the gene expression of two cohorts treated with immunotherapy. The first cohorts (GSE35640) included patients with metastatic melanoma treated with MAGE-A3 immunotherapeutic (18). Clinical benefit included objective responders (complete and partial) according to RECIST 1.0 (19). Only 22 responders and 34 nonresponders were used in our analysis. Another cohort (GSE78220) included melanoma biopsies treated with anti–PD-1 therapy (20). We analyzed gene expression of 5 complete responders and 13 progressive disease samples.

All cell lines used in the study were obtained from ATCC and were authenticated by short tandem repeat (STR) profiling and tested for the absence of Mycoplasma. Cells were passaged twice from collection or thawing for use in experiments. In vitro cell viability assays for mutations were adapted from one of our recent studies (21). Two cell lines Ba/F3 and MCF10A were used and mutation functional calls, including activating, inactivating, inhibitory, noninhibitory, and neutral, were determined by comparing with the corresponding wild-type counterparts. We used the consensus calls that were consistent in both cell lines.

ESTIMATE (Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data) was used to estimate the immune score, stromal score, and tumor purity of each patients across 33 cancer types (22). The FPKM-based gene expression profiles were used as input of this method. In addition, TIMER (https://cistrome.shinyapps.io/timer/) was used to estimate the abundances of member cell types in a mixed cell population, using gene expression data (23).

Pathway mutation burden (PMB) score was defined to evaluate whether the immune pathways were likely to mutate in a specific cancer type. This score considers the proportion of mutated genes as well as the coverage of patients in each cancer type. For an immune pathway i in cancer type j, the score is defined as following:

where |{m_{ij}}$| is the number of mutated immune genes in pathway i; |{m_i}$| is the total number of genes in pathway i; |{n_{ij}}$| is the number of patients of cancer j that with gene mutations in pathway i; and |{n_j}$| is the total number of patients sequenced in cancer j.

Here, we considered three immune response–related scores estimated from gene expression. The immune score that represented the infiltration of immune cells in tumor tissue was assessed based on ESTIMATE (22). The MHC score was formulated from the gene expression levels of the “core” MHC-I set (including HLA-A, HLA-B, HLA-C, TAP1, TAP2, NLRC5, PSMB9, PSMB8, and B2M; ref. 24). The FPKMs of genes were first log-transformed and then median-centered. The mean expression levels of these core genes were defined as the MHC score. The cytolytic activity (CYT score) was based on transcript levels of two key cytolytic effectors, granzyme A (GZMA) and perforin (PRF1), found in previous studies (25).

To investigate whether the immune-related scores could be used to classify the patients with different immune responses to immunotherapy, we obtained gene expression profiles of patients with metastatic melanoma treated with MAGE-A3 immunotherapy (25) or anti–PD-1 (20). The three types of immune response–related scores were calculated for each patient. Next, we used the ROC curve analysis and precision–recall curve to evaluate the effects of these scores. This process was performed by using “ROCR” package in R program. Additional analysis was performed using patients from the TIDE server (26).

The domainXplorer algorithm was used to identify the protein regions that are associated with immune scores (27). The protein functional regions were defined as protein domains in Pfam or intrinsically disordered regions identified by Foldindex. In addition, potential domains identified by AIDA (28) were also included in our analysis. Here, we only focused on immune genes–related functional regions.

In brief, domainXplorer used three statistical tests to evaluate whether mutations in a functional protein region were correlated with the immune response–related scores (27). The first linear model was defined as follows:

where “ImS” is the immune score of each patient, “T” is the tissue of origin of each patient (here is the cancer type), and “D” indicates whether the functional region is mutated or not. In the next step, Wilcoxon rank sum test was used to compare the immune response–related scores in patients with mutations in the functional region being analyzed against those with mutations in other regions of the same immune gene (Test-2) or those without mutations in the immune gene at all (Test-3). Protein regions with p1 < 0.05, p2 < 0.05, and p3 < 0.05 were identified as to be associated with immune response. In addition, we also identified the protein regions that were correlated with the proportion of immune cell infiltration based on the same method.

To identify differentially expressed genes in each cancer type, we downloaded raw RNA-seq counts for all tumor samples, compared with normal tissues. We then used Limma-Voom (29, 30) to identify the genes and considered the tumor purity as a factor. Genes with P values less than 0.05 and 4-fold changes in expression were identified as differentially expressed genes in each cancer type.

To evaluate the similarity of cancer types based on the differentially expressed immune genes, we calculated the Jaccard index for each pair of cancer

where |DE{G_i}$| and |DE{G_j}\ $|are the differentially expressed immune genes in cancer type i and j. The directions of differentially expressed genes were not considered in our analysis.

Gene set enrichment analysis (GSEA) was used to identify the enriched immune pathways by differentially expressed genes in each cancer type (31). First, genes were ranked by negative log₁₀ of the differential expression analysis–derived FDR multiplied by the sign of the logFC (log fold change). The “weighted” enrichment statistics were calculated for enrichment or depletion of each immune pathway in specific cancer type.

The somatic mutations were collected from TCGA project and we obtained the HLA alleles for tumor samples from TCIA (https://tcia.at/home; ref. 32). Next, NetMHCpan 4.0 was employed for neoantigen prediction based on the information on mutations and HLA alleles (33). All the wild-type and mutated peptides with 8–11 amino acids were extracted and predicted for the binding affinity. Rank of the predicted affinity was compared to a set of random natural peptides. Rank threshold for strong binding peptides was 0.5% and for weak binding peptides was 2.0%. A cutoff of 0.5% was used to define mutations giving rise to neoantigens. The predicted binding affinity was in nanomolar units.

Human protein–protein network was integrated into our analysis to prioritize immune dysregulated genes. First, we extracted the subnetwork formed by immune-related genes from an integrated functional network (34). The main component of this subnetwork consists 6,060 interactions among 1,202 immune-related genes. Next, the immune genes were ranked based on the correlation between gene expression and immune scores, MHC scores, and CYT scores. Top 50 genes with higher positive or negative correlation with immune, MHC, and CYT scores were selected as seed genes. In addition, the genes in which mutations are associated with immune-related scores were also integrated as seed genes. To simulate the propagation of the genetic alterations through the network, we employed the random walk with restart (35). When an immune gene in the network is ranked in top 50, a value “1” is initially assigned to the gene. Then the genetic alterations are propagated along the neighbors such that the genetic alteration probability was calculated for all genes in the network according to the following equation:

where |{P_0}$| is the initial binary probability vector, A denotes a degree-normalized adjacency matrix of the immune subnetwork, and |\alpha $| determines the diffusion degree of the genetic influence through the network. We used the optimal value (⁠|\alpha = 0.7$|⁠) in our analysis (36). After numbers of iteration, we obtained the final probability for each gene in the network. Genes with greater score than average probability were identified as immune dysregulated genes for each cancer type.

Next, we integrated drug-associated gene signatures to identify the drug candidates for each cancer type. We first obtained the expression signature genes from The Connectivity Map (also known as CMAP; ref. 37). For each drug, there are two types of signature genes, one is the upregulated genes and the other one is the downregulated genes. We used hypergeometric test to evaluate whether the drug-perturbed gene signatures were enriched in cancer-specific immune gene sets. For each cancer type, we obtained the P value and then all the P values across 33 cancer types were combined to a combined P value.

where w is the weight for the individual P values, and z is the Z-score of P values. Here, we set the same weights for all individual P values. The drugs with adjusted P values less than 0.01 were identified. We obtained the drugs that had upregulated/downregulated gene signatures enriched for genes positively/negatively correlated with immune scores.

The patients were divided into two groups based on the median expression of specific immune gene. Log-rank test was used to evaluate the survival difference between two groups. Genes with P values less than 0.05 were considered as statistically significant.

TMB, a measurement of the overall number of genetic alterations observed in a cancer sample, has been suggested as a potentially prognostic marker for immunotherapy (2). In addition, several gene signatures have also been proposed to predict the immune response of patients (Fig. 1A). We first evaluated whether three distinct proposed immune gene expression signatures were associated with immunotherapy response: (i) immune score, which is a measure of total immune infiltrate into the tumor; (ii) MHC score, which is a measure of antigen presentation required for tumor cell recognition by T cells and subsequent T cell–mediated killing; and (iii) cytolytic activity, which is a measure of cytolytic enzymes used by immune cells to kill tumor cells. Recombinant MAGE-A3 protein combined with immunostimulant AS02B or AS15 has been administered to patients with early metastatic melanoma (18). On the basis of the gene expression data in response to treatment with MAGE-A3 immunotherapy, we found that the patients that responded to immunotherapy had significantly higher immune scores (P = 0.003), MHC scores (P = 1.7e-4), and cytolytic activity (P = 0.007, Fig. 1B). Moreover, we found that these three immune-related scores indeed distinguished the responder from non-responder patients with area under the curve of the ROC (AUROC) from 0.71 to 0.79 (Fig. 1C). To further validate this, we analyzed another gene expression dataset for patients with melanoma treated with anti–PD-1 immune checkpoint blockade (20). We also found that patients that were responsive to the PD-1 immunotherapy exhibited higher immune-related scores (Supplementary Fig. S1A). The AUROCs for classification ranged from 0.52 to 0.85, and the precision–recall curves showed clear separation of responders versus nonresponders (Supplementary Fig. S1B–S1D), which was similar to the results obtained from gene signatures in previous studies (38). These results suggest that the immune-related scores can potentially reflect the immune milieu of tumors in patients and help identify specific factors that may modulate these environments.

We next explored the generalizability of these scoring systems in predicting patient immune responses across cancer types. We first calculated the three immune-related scores (immune score, MHC score, and cytolytic activity) for all tumor patients from the TCGA project (Fig. 1D). We found that the distributions of the immune-related scores were diverse across 33 cancer types. For example, tumors from patients with kidney renal clear cell carcinoma (KIRC) exhibited higher immune-related scores, and immune checkpoint blockade has shown promising results for first-line treatment of KIRC in multiple phase III clinical trials (39). Moreover, although the three immune-related scores appeared to be correlated with each other, the Pearson correlation coefficients ranged from 0.51 to 0.79 (Fig. 1E). This observation implied that these scores might be complementary with each other; therefore, integrating these scores would facilitate the identification of novel candidate gene targets for immunotherapy.

It has been suggested that the mutations in immune pathways could impact tumor-immune interactions (40). It is unclear, however, how the immune mutational burden (IMB) is associated with patient outcome. Moreover, it seems that the mutation burdens across pathways are varied (11). Therefore, it remains elusive if and how mutations in immune-related pathways contribute to higher mutation burden. To systematically investigate the global functional genomic landscape of immune-related genes, we manually curated 2,273 immune-related genes from literature (12, 13). These genes were further classified into 18 immune-related pathways (Fig. 2A). Next, we analyzed whole exomes from 100 melanoma patients treated with CTLA-4 blockade (41). We first calculated the number of mutations (TMB) and the number of mutations located in immune-related genes (defined as IMB). We further compared TBM and IMB between responders and nonresponders. The mutation burden and immune burden of responders were significantly higher than those of nonresponders (Supplementary Fig. S2A; P < 0.05). In addition, mutation burden and immune burden could predict immune response with similar power (Supplementary Fig. S2B; AUCs = 0.68). Furthermore, we calculated the IMB for each patient in TCGA cohort. We found that the cancer types that are likely to respond well to immunotherapy had higher proportion of immune-related gene mutations, such as skin cutaneous melanoma (SKCM) and kidney cancer (Supplementary Fig. S2C). These results suggest that the immune burden is a useful indicator of immune response; target sequencing of immune-related genes would help predict the immune response of patients.

Next, we defined a PMB score to evaluate to what extent each immune-related pathway is perturbed by genomic alterations in tumor. This score is calculated based on the mutation frequency and the proportion of mutated genes in each pathway. We found that the essential genes identified by CRISPR-Cas9 KO screens (13) had higher PMB scores across cancer types (Fig. 2B). Moreover, we observed another pathway-cytokine receptors (CR) with higher PMB scores in LUSC, SKCM, and DLBC (Fig. 2B). We next focused on the CR mutations in LUSC (Fig. 2C). Consistent with previous studies, the mutations in CR pathway exhibited a mutually exclusive pattern (42). Moreover, the mutations in CR pathway had higher functional impact scores than other mutations evaluated by metaSVM (17) (Supplementary Fig. S2D). The top-1 gene (PLXNA4) ranked by mutation frequency was mutated in approximately 19% of patients with lung cancer. PLXNA4 was shown previously to promote tumor progression and tumor angiogenesis by activating VEGF and FGF signaling (43). We further found that the expression of PLXNA4 was correlated with the patients' survival in LUSC (Supplementary Fig. S2E; P = 0.017). These results suggest that the mutations in cytokine receptors pathway might play critical roles in cancer and the proposed PMB index could help prioritizing key mutations in cancer pathways.

Next, we interrogated if and to what extent the human immunome is altered in cancer at the gene expression level. We identified differentially expressed genes across 18 cancer types with more than five normal samples. As tumor purity may affect the identification of differentially expressed genes (44), we first estimated the tumor purity of each sample (Supplementary Fig. S3A) and identified the differentially expressed genes using tumor purity as a covariable in the regression model. We found that the immune-related genes were significantly perturbed in cancer when compared to other coding genes (Fig. 3A; Supplementary Fig. S3B; P < 0.05 for all cancer types, Fisher exact test). The odds ratios (OR) of immune genes versus nonimmune genes were higher than 1.0 in all cancer types. These results indicate that these genes could play a role in regulating tumor growth.

Cancer types with similar tissue origins often show similar transcriptome profiles (45). Next, we calculated the Jaccard index for each pair of cancer based on the similarity of differentially expressed immune-related genes (see details in Materials and Methods). We found that cancer lineages with similar tissue origin were clustered together (Fig. 3B), such as lung adenocarcinoma (LUAD) and lung squamous cell carcinoma (LUSC), cholangiocarcinoma (CHOL) and liver hepatocellular carcinoma (LIHC), and three types of kidney cancer. These results indicate that these cancer types show immunome profiles within similar tumor microenvironment.

In addition, we ranked all genes based on their differential expression and used GSEA to explore which immune-related pathways were likely to show transcriptional perturbation. Sixteen of 18 immune-related pathways were significantly altered in at least one cancer type (Fig. 3C). The antigen processing and presentation, cytokines, and cytokines receptors pathways were perturbed in more than 16 cancer types (Supplementary Table S1). Moreover, we found that these immune-related pathways showed consistent trend (enriched or depleted) across cancer types. For instance, although tumors have been observed to downregulate antigen presentation for immune escape (46), we found antigen processing and presentation pathway was likely to be activated in breast-invasive carcinoma and glioblastoma (BRCA and GBM; Fig. 3C and D; Supplementary Fig. S3C). In addition, cytokines and cytokine receptors pathways were repressed in BRCA and GBM (Fig. 3E and F; Supplementary Fig. S3D; FDR < 0.001). Moreover, we performed expression analysis based on immunologic gene signatures in MSigDB (47), which provides more specific immune pathways. In total, 4,872 immune-related gene signatures were analyzed. We found that approximately 67.3% of these immunologic gene signatures were enriched in differentially expressed genes in at least one cancer type (Supplementary Fig. S4). Moreover, there were 27 gene signatures significantly enriched in >11 cancer types (Supplementary Table S2). These results suggest that immunologically relevant genes are likely to be perturbed in cancer.

Growing evidence suggests that mutations in certain genes influence the immune response (27). Identification of these genes may provide valuable insights for improving immunotherapy. Thus, we explored “domainXplorer” to identify critical regions of immune genes that are associated with immune response in cancer (27). Given the importance of the proposed immune-related scores in predicting immune response, we next determined whether cancer mutations in candidate driver regions were significantly correlated with these scores (Supplementary Table S3). Our analysis yielded a total of 209 protein regions in 145 immune-related genes that are likely associated with cancer immunity (Fig. 4A). Specifically, mutations in seven genes (UBXN1, LTBP2, STAT1, NGFR, NRAS, MET, and LTBP4) were significantly correlated with all of three types of immune-related scores. LTBP4, a member of the latent TGFβ binding protein gene family, had been linked to several human diseases (48). We found that patients with EGF_CA domain mutations exhibited significantly higher immune scores, MHC scores, and cytolytic activity (Fig. 4B, P < 0.05). Moreover, we also identified mutations of STAT1 significantly correlated with these scores, which was consistent with previous observation that STAT1 played an important role in the innate immune response (49). These results collectively expand the catalog of potential cancer immune drivers and highlight the importance of taking into account the protein structural context to identify patients that are likely to clinically benefit from immunotherapy.

Our above analyses identified candidate protein regions that are potentially correlated with immune response. However, we are still lack of knowledge about the potential mechanism. An emerging role in treatment response has been attributed to immune cell infiltration in human tumors (7). Thus, we next explored to what extent the mutations in these regions were correlated with immune cell infiltration. We first identified the protein regions whose mutations were significantly correlated with the abundance of six tumor-infiltrating immune cells (TIIC) subsets (B cells, CD4 T cells, CD8 T cells, macrophages, neutrophils, and dendritic cells; ref. 50). We also identified the proteins that were associated with B cell receptor (BCR) and TCR diversity, by CPK (CDR3s per kb of reads) that evaluates the clonotype diversity (51). We found that both MHC and CYT scores showed significant association with the expansion of tumor-infiltrating B cells after correction for tumor purity (Supplementary Fig. S5A). In addition, MHC score was weakly associated with TCR diversity, potentially due to clonal expansion of tumor antigen-specific T cell populations. CYT score on the other hand was strongly correlated with CPK (Supplementary Fig. S5B). We found that the protein regions that were correlated with immune score, MHC score, and cytolytic activity were significantly overlapped with those correlated with TIIC abundance (Fig. 4C). Especially, there were a higher number of protein regions correlated with CD8 T cells, neutrophils, and dendritic cells.

Specifically, we found that seven genes, whose mutations correlated with all three types of immune-related scores, were associated with at least one type of immune cell abundance (Fig. 4D). For example, patients with LTBP4 mutations exhibited significantly higher abundance of CD8 T cells and neutrophils (Fig. 4E). CD8⁺ T cells are crucial mediators of antitumor immune response and the targets of checkpoint blockade (52), while the role of neutrophils is complicated as they can exert either tumor enhancing or tumor clearance effect (53). Moreover, we found that not only the genomic variants in LTBP4 were correlated with immune cell abundance, but patients with copy number alterations also exhibited significantly different abundance (Supplementary Fig. S5C). These results suggest that the genomic alterations might contribute to immunotherapy response through affecting the immune cell infiltration in cancer.

Increasing studies have demonstrated that mutation-derived neoantigens form a major “active ingredient” of successful cancer immunotherapies (3). Next, we predicted the neoantigens that were derived from mutations of TCGA project based on NetMHC (33). There was a limited proportion of neoantigens resulting from mutations in protein coding genes and approximately 1.2% of all mutations in immune-related genes could potentially generate neoantigens (Fig. 5A). This is consistent with a recent study showing that neoantigens arose from mutations in only about 60% of the genes (54), while the majority of neoantigens might originate from alternative or noncoding sources. We found that the immunologically relevant genes (from ImmPort) tended to generate significantly more neoantigens (such as PIK3CA, MUC4, EGFR) compared with coding genes in general (Fig. 5A, P = 1.18e-4, Fisher exact test). Approximately 25% of the mutations in PIK3CA and HRAS could generate neoantigens (Fig. 5B). This result is consistent with the observation that RAS plays a role in immunotherapy (55).

Next, we focused on mutations in PIK3CA and EGFR that give rise to neoantigens. On the basis of the viability-based functional screen in two cell lines (Ba/F3 and MCF10A; ref. 21), we found that the mutations that encode neoantigens were likely to be activating mutations in both cell lines (Fig. 5C and D). For example, we identified 46 missense mutations that could generate neoantigens in PIK3CA. Among these mutations, 41 mutations were in the functional dataset. We found that 40 mutations were activating mutations and one was inactivating mutations. Moreover, we identified a mutation cluster (including R38C/H, E81K, R88Q, R93W, K111E, and R115L) in PIK3CA (Fig. 5E). For EGFR, we identified 17 mutations that could generate neoantigens. In total, 14 mutations were screen for function and 11 mutations were activating, two were neutral and 1 inconclusive in two cell lines (Fig. 5D). We also found seven mutations formed a cluster in 3D (including R108K, R252C/P, and A289V/D; Fig. 5E). We hypothesize that functional mutations that cause a conformational change in the protein may be more likely to be recognized by MHC molecules for immune pruning, but if these mutations subsequently suppress inflammatory signaling, it may allow for proliferation of tumor cells with activating mutations.

Our above results indicate widespread genomic and transcriptome perturbations in the cancer immunome network. An integrated landscape of these perturbations offers a basis for further advancing our understanding of the immune regulation in cancer. Thus, we proposed a Network-based Integrative model to Prioritize Potential immune respondER genes (NIPPER Fig. 6A). Briefly, we first identified the genes whose mutations or expression were significantly correlated with immune scores, MHC scores, and CYT scores. These genes were then mapped to an immune-related protein–protein interaction (PPI) network, which consisted of 6,060 interactions among 1,202 immune-related genes. Next, the mutational effects were propagated in the PPI network. Finally, the genes with driver probability greater than expected were identified as immune response-related signature genes. This process was repeated for all cancer types. Here, the genes positively/negatively correlated with immune-related scores were analyzed, separately.

On the basis of the proposed network-based model, we identified 39 genes as signature genes in 33 cancer types (Fig. 6B). High expressions of these genes were likely to be correlated with immune-related scores. These genes were involved in 11 immune-related pathways, and 8 genes (such as STAT1, B2M, TAP1, and TAP2) had been identified as essential genes in immune response. Moreover, we found that these candidate immune response-related genes were clustered together as modules in the PPI network (Fig. 6C). Seven of the 8 essential immune genes were linked in a dense module, and STAT1 was clustered together with more novel immune response-related genes. Interestingly, we identified one candidate gene-PDIA3, which linked 7 known essential genes. The expression of this gene was correlated with the patients' survival in several cancer types, including CESC, LUAD, HNSC, and GBM (Fig. 6D; Supplementary Fig. S6). These observations suggest that this gene may play critical roles in immunotherapy, which is consistent with the conclusion of one recent study (56).

Next, we validated these genes in a cohort of melanoma patients that underwent treatment with adoptive T-cell transfer (24). We found that the expression of these genes could classify the responder and nonresponder patients with an AUC of 0.713 (Fig. 6E). This is similar to the results based on PD-1 or PD-L1 expression. In another cohort of patients with melanoma treated with anti–PD-1, we found that the power of NIPPER was higher than PD-1 and PD-L1 expression-based and similar as TMB-based methods (Supplementary Fig. S7). We also used TIDE to compare our gene signatures with other public signatures in 17 datasets (26). Our NIPPER signatures achieved AUCs > 0.6 in most (9/15) datasets (Supplementary Fig. S8). Moreover, the average expression of these genes was significantly associated with improved overall survival (Fig. 6F, log-rank P = 0.02). Taken together, these results indicate that signaling networks identified by NIPPER are relevant to response to multiple immunotherapy modalities, and may offer insight on critical signaling cascades that can be targeted to improve immunotherapy outcomes.

Next, we evaluated mutations in the genes identified by NIPPER in different cancers. Here, all mutations in prioritized genes across cancer types were analyzed. We found that the mutations in these genes were likely to alter the binding affinity of HLAs (Fig. 7A, P < 0.001, Wilcoxon rank sum test). Specifically, we identified 135 mutations in 28 genes that perturbed binding of HLAs (Fig. 7B; Supplementary Table S4). These mutations were distributed in various types of cancer. We identified six mutations in PIK3R1 (encoding p85α) that reduced binding of HLAs (Fig. 7B) and some mutations clustered in a possibly functional spot (Fig. 7C). Consistently, we identified two mutations (N564D and K567E) that increased growth of cancer cell lines (Fig. 7D). These results suggest that our network-based method helps identify cancer functional mutations by their change in binding affinity to HLAs.

Identification of candidate small-molecule drugs is critical for the immunotherapy in cancer. Herein, we identified candidate drugs that could increase the expression of identified immune gene signatures based on GSEA analysis (Fig. 6A). We obtained more than 6,000 small molecules and their perturbation gene signatures from The Connectivity Map (57). We identified 49 small molecules that significantly perturbed the expression of immune gene signatures (P < 0.05, Fig. 7E; Supplementary Table S5). Interestingly, we identified propofol as the top one small molecule for increased immune response. Evidence had shown that propofol could possibly induce a favorable immune response in terms of preservation of IL2/IL4 and CD4⁺/CD8⁺ T-cell ratio in the perioperative period for breast cancer (58). In vitro and in vivo studies had also suggested that propofol could independently reduce migration of cancer cells and metastasis (59). Taken together, these results suggest that the network-based integrative model not only can identify immune response–related gene signatures, it also helps prioritize candidate drugs for immunotherapy in cancer.

Systematically understanding the regulation of immune systems opens more avenues for cancer immunotherapy with a potent clinical efficacy. By integration of pan-cancer omics datasets, we have examined human tumor correlations in the context of immune-related pathways to an extent that is not previously possible. High-throughput gene expression data have been widely used to investigate differentially expressed genes in various types of cancer. Strikingly, we observed extensive immunologic gene signatures with widespread perturbation in cancer compared with normal samples. These deregulated immunologically relevant genes are enriched in “antigen processing and presentation,” “cytokines,” and “cytokines receptors” pathways. Cytokines are critical molecular messengers for cells of the immune system to communicate with one another, which can directly stimulate immune effector cells and enhance tumor cell recognition. Increasing studies have demonstrated that cytokines have broad antitumor activity (60). Our analysis revealed that there are widespread perturbations of cytokine pathways in more than 80% of cancer types. A more detailed understanding of cytokine pathway regulation will provide new opportunities for improving cancer immunotherapy.

It is still unclear which factors are key players in regulating the immunome homeostasis. Thus, we systematically analyzed two types of genetic regulation (including CNVs and eQTLs) in immune genes (Supplementary Fig. S9A–S9C). We found that a number of immunologically relevant genes recurrently dysregulated in various types of cancer. For example, upregulation of BIRC5 is correlated with CNV amplification (61), and downregulation of CCL14 is correlated with CNV deletion across cancer types. Moreover, eQTLs analysis revealed a hotspot locus that modulates the expression of DEFB1 (Supplementary Fig. S10A–S10D), which might function through perturbation of CTCF binding. Although a number of immune-related genes show correlation between copy number variations/mutations and gene expression, the relationship is weak and could not fully explain immune-related transcriptome perturbations. Thus, it is likely that other genetic (such as posttranscriptional regulation) or epigenetic variables (such as DNA methylation or histone modification) also contribute to the transcriptome perturbations in cancer.

Prior studies have proposed that the immune score, MHC score, and CYT score are correlated with immune response. For example, Rooney and colleagues quantified the cytolytic activity and identified associated properties across 18 tumor types (25). We found that these three scores that interrogate different aspects of immune signaling are complementary with each other, providing more information than any one signature alone. Further information was gained by utilization of a network-based approach, specifically a random walk with restart network propagation, which considers both the mutations in individual genes as well as the topology of interactions between the corresponding proteins. Integration of the network topology information has been demonstrated to have much higher accuracy in identifying cancer-related genes (62). Thus, we proposed an integrated network model to predict the critical signaling cascades for immunotherapy response. This method identified several well-known immune response-related genes, such as B2M (63), TAP1 (64), and TAPBP (65), as well as several novel genes such as PDIA3. Moreover, we predicted candidate small molecules that may sensitize to immunotherapy based on whether the treatment could induce the expression of these signature immune genes. In particular, propofol was identified as one of the top-ranked drug candidates for immunotherapy. Propofol is most commonly used in clinical anesthesia and propofol exhibits a good inhibitory effect on tumor recurrence and metastasis (59). Our current study found that propofol likely induces the expression of immune response-related genes, which may contribute to improved efficacy of immunotherapy.

Taken together, our integrative analysis of the human immunome across cancer types reveals several candidate gene and mutation markers that may be critical for understanding patient responses to immunotherapy. Future studies will need to evaluate the potential molecular functions of our identified signaling cascades by low-throughput experiments.

D.J. McGrail reports grants from NCI (NCI grant K99CA240689) during the conduct of the study. S.A. Shukla reports grants from NCI during the conduct of the study, nonfinancial support and other compensation from Bristol-Myers Squibb (equity), personal fees from Neon Therapeutics (consulting), other compensation from Agenus Inc. (equity), Agios Pharmaceuticals (equity), BreakBio Corp (equity), and Lumos Pharma (equity) outside the submitted work. C.J. Wu reports other from Neon Therapeutics (founder and SAB member) outside the submitted work, as well as a patent for PCT/US2011/036665 licensed to Neon Therapeutics. G.B. Mills reports personal fees, nonfinancial support, and other compensation from AstraZeneca (clinical trials support), personal fees and nonfinancial support from Chrysalis, personal fees, nonfinancial support, and other compensation from ImmunoMet (stock), personal fees and nonfinancial support from Ionis, personal fees, nonfinancial support, and other compensation from Lilly (clinical trials support), personal fees and nonfinancial support from PDX Pharma, personal fees, nonfinancial support, and other compensation from Signalchem (stock), personal fees and nonfinancial support from Symphogen, personal fees, nonfinancial support, and other compensation from Tarveda (stock), personal fees and nonfinancial support from Turbine and Zentalis, other compensation from Catena (stock), Genentech (clinical trials support), GSK (clinical trials support), grants from Adelson Medical Research Foundation, Breast Cancer Research Foundation, Komen Research Foundation, Ovarian Cancer Research Foundation, and Prospect Creek Foundation during the conduct of the study. N. Sahni reports grants from Breast Cancer Alliance, Ovarian Cancer Research Alliance, Alfred P. Sloan Foundation, and Cancer Prevention and Research Institute of Texas during the conduct of the study. S. Yi reports grants from Susan Komen Foundation and grants from NCI during the conduct of the study. No potential conflicts of interest were disclosed by the other authors.

Y. Li: Conceptualization, data curation, validation, investigation, visualization, methodology, writing-original draft. B. Burgman: Validation, investigation, methodology, writing-review and editing. D.J. McGrail: Conceptualization, funding acquisition, validation, investigation, methodology, writing-original draft, project administration, writing-review and editing. M. Sun: Conceptualization, resources, supervision, funding acquisition, validation, investigation, methodology, writing-original draft, project administration, writing-review and editing. D. Qi: Validation, investigation, methodology, writing-review and editing. S.A. Shukla: Resources, data curation, validation, investigation, methodology, writing-review and editing. E. Wu: Writing-review and editing. A. Capasso: Project administration, writing-review and editing. S.-Y. Lin: Resources, data curation, investigation, writing-review and editing. C.J. Wu: Resources, data curation, project administration, writing-review and editing. S.G. Eckhardt: Funding acquisition, project administration, writing-review and editing. G.B. Mills: Resources, data curation, funding acquisition, validation, methodology, writing-review and editing. B. Li: Resources, funding acquisition, validation, methodology, project administration, writing-review and editing. N. Sahni: Conceptualization, funding acquisition, validation, methodology, writing-original draft, project administration, writing-review and editing. S.S. Yi: Conceptualization, resources, supervision, funding acquisition, validation, investigation, methodology, writing-original draft, project administration, writing-review and editing.

The authors would like to acknowledge Kevin Drew and Edward M. Marcotte from University of Texas at Austin for helpful discussions. The authors are grateful to contributions from TCGA Research Network Analysis Working Group, and also acknowledge the Biomedical Research Computing Facility at UT Austin and Texas Advanced Computing Center (TACC) for high-performance computing assistance. This work was supported by NIH grant (K22CA214765 to S.S. Yi), Komen Foundation grant (CCR19609287 to S.S. Yi), Young Investigator Grant by Breast Cancer Alliance (to N. Sahni), Liz Tilberis Early Career Award by Ovarian Cancer Research Alliance (Grant# 649968 to N. Sahni), Alfred P. Sloan Scholar Research Fellowship (FG-2018–10723 to N. Sahni), Cancer Prevention and Research Institute of Texas grants RR170079 (to B. Li) and RR160093 (to S.G. Eckhardt). N. Sahni is a CPRIT Scholar in Cancer Research with funding from the Cancer Prevention and Research Institute of Texas (CPRIT) New Investigator Grant RR160021. D.J. McGrail is supported by NCI grant K99CA240689. S.A. Shukla acknowledges support from the NCI (R50RCA211482). A. Capasso is supported by Department of Defense grant CA191245. C.J. Wu is a Scholar of the Leukemia and Lymphoma Society. G.B. Mills is supported by a kind gift from the Miriam and Sheldon Adelson Medical Research Foundation, SAC110052 from the Susan G Komen Research Foundation, Breast Cancer Research Foundation, Ovarian Cancer Research Alliance, Prospect Creek Foundation, and NCI grant (U01CA217842).

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