Multi-omics analysis reveals neoantigen-independent immune cell infiltration in copy-number driven cancers

Abstract

To realize the full potential of immunotherapy, it is critical to understand the drivers of tumor infiltration by immune cells. Previous studies have linked immune infiltration with tumor neoantigen levels, but the broad applicability of this concept remains unknown. Here, we find that while this observation is true across cancers characterized by recurrent mutations, it does not hold for cancers driven by recurrent copy number alterations, such as breast and pancreatic tumors. To understand immune invasion in these cancers, we developed an integrative multi-omics framework, identifying the DNA damage response protein ATM as a driver of cytokine production leading to increased immune infiltration. This prediction was validated in numerous orthogonal datasets, as well as experimentally in vitro and in vivo by cytokine release and immune cell migration. These findings demonstrate diverse drivers of immune cell infiltration across cancer lineages and may facilitate the clinical adaption of immunotherapies across diverse malignancies.

Fig. 1

Integrative multi-omics network analysis framework to reveal determinants of immune infiltration across human cancers. a Flow chart for analysis of pan-cancer indicators of immune invasion. Cytotoxic T-cell (CTL) levels were determined from proteomic data and integrated with gene expression, genetic mutations, phospho-/total- proteomics, and interactome networks. Identified markers were then expanded pan-cancer, and potential novel drivers were experimentally verified. b Protein data sets and classification of CTL high patients. CTL high breast cancer patients were classified by positive expression of the cytolytic enzyme granzyme B and CD8 CTL surface marker, and also showed enrichment of the cytolytic enzyme perforin. P values were determined by Wilcoxon rank-sum test. Box indicates median with interquartile range, and whisker length determined by the Tukey’s method

Fig. 2

Proteomic identification of patients with high CTL invasion reveals that breast cancer invasion does not correlate with neoantigens. a Receiver-operator characteristic plot for prediction of CTL high patients based on neoantigen levels in colorectal cancer patients. Area under the curve (AUC) equal to 0.757 (N = 64). b Receiver-operator characteristic plot for prediction of CTL high patients based on neoantigen levels in breast cancer patients. Area under the curve (AUC) equal to 0.532 (N = 80). c CTL score calculated across 19 cancers, plotted in order of median for each cancer. Dots represent individual patients, and lines are median with interquartile range. See Supplementary Fig. 3. Sample sizes given in Supplementary Table 1. d Neoantigens per patient across 19 cancers, plotted in order of median CTL score for each cancer. Dots represent individual patients, and lines are median with interquartile range. Sample sizes given in Supplementary Table 1. e Plot of CTL score as a function of neoantigen load across all patients (N = 7835). Inset r is Spearman’s correlation coefficient and corresponding P value. f Pan-cancer analysis of neoantigen-CTL Spearman correlation, with size of dots topping each bar representing significance level (P). Sample sizes and exact P-values given in Supplementary Table 2

Fig. 3

Identification of cytokines responsible for CTL invasion across multiple cancers. a Gene set enrichment analysis for total proteins in CTL high breast cancer. b Leading edge analysis of top ten enriched gene sets shows soluble factors are driving most of the enriched gene sets. c Volcano plot showing differential expression of secreted proteins in CTL high vs. CTL low patients. Dotted line indicates false discovery rate (FDR) of 5%. Proteins which show positive correlation between protein and gene transcript levels are indicated in blue. See Supplementary Fig. 4. d Pan-cancer analysis of correlation between CTL score and cytokine expression levels. Correlation coefficient (R) indicated by color of the spot, with significant correlations indicated by light green background. Sample sizes and exact P-values given in Supplementary Table 3

Fig. 4

ATM as a driver of CTL invasion in breast cancer. a Gene set enrichment analysis for differentially expressed phospho-proteins shows enrichment of the ATM pathway. b Gene set enrichment analysis for breast cancer CTL high gene signature (see Supplementary Fig. 5) shows enrichment of ATM related genes. c Receiver-operator characteristic (ROC) plot for prediction of CTL high patients based on ATM phospho-S1981 levels in breast cancer patients. Area under the curve (AUC) equal to 0.82 (N = 80). d Correlation of total ATM with ATM phosphorylated at serine 1981 (N = 80). Spearman correlation coefficient (r) and corresponding P-value indicated above plot. e ROC plot for prediction of CTL high patients by total ATM protein (N = 80). f Correlation of CTL Score with total ATM determined by RPPA for the TCGA breast cancer cohort (N = 892). Spearman correlation coefficient (r) and corresponding P-value is indicated in the above plot. g Analysis of breast cancer patients with genomic alteration of ATM shows loss of correlation between ATM and CTL invasion compared to patients with wild type (WT) ATM. Correlation coefficient distributions were determined by random subsampling of each population to ensure equal sample size. Error bar represents 95% confidence interval, P-values determined by Wilcoxon rank-sum test

Fig. 5

Enrichment of transcription factors for CTL cytokines in the ATM network. a Work flow to identify transcription factors within the ATM second neighbor network. b Enrichment of transcription factors that drive cytokines implicated in CTL recruitment. Blue arrows indicate number of transcription factors in ATM second neighbor network. The null distribution shown in the background represents number of transcription factors bound in randomly selected gene sets of equal size, which was used to determine an empirical P-value (inset). c Correlation of ATM phosphorylated at S1981 and cytokine protein levels in PDX breast cancer models (N = 27). Spearman correlation coefficient (r) and corresponding P-value indicated above plot. Error bars represent mean ± SEM. d Analysis of breast cancer patients with genomic alteration of ATM shows loss of correlation between ATM and cytokine expression compared to patients with wild type (WT) ATM. Correlation coefficient distributions were determined by random subsampling of each population to ensure equal sample size. Error bar represents 95% confidence interval, P-values determined by Wilcoxon rank-sum test. Error bars represent mean ± SEM

Fig. 6

Phosphorylation of ATM induces cytokine secretion and PBMC migration. a Quantification of CCL5, CXCL10, and IL16 secretion by ELISA following irradiation of BT549 breast cancer cells either in presence or absence of ATM inhibitor KU-55933 at 10 μM. Values plotted as mean ± std of duplicate spots. See Supplementary Fig. 7. b Schematic for collection of conditioned media for PBMC transwells. Final concentration of ATM inhibitor KU-55933 was equal in all conditions (10 μM). c Migration of PBMCs towards tumor cell conditioned media from triple negative breast cancer cell lines MDA-MB-231 and BT-549 and luminal breast cancer cell lines ZR-75-1. Error bars represent ± standard deviation of duplicate runs. Significance determined by ANOVA using a Holm-Sidak post-hoc test

Fig. 7

ATM phosphorylation correlates with CTL invasion and cytokine expression in a murine breast cancer model. a Schematic of murine derived syngeneic transplant (MDST) model. b Representative CD8 (green) immunostaining with DAPI nuclear counterstain (blue) of MDST tumor sections for a CTL-high MDST (i) and CTL-low MDST (ii). Scale bar = 50 μm. c Correlation of ATM pS1981 determined by RPPA with CD8 immunostaining (N = 11 MDST models). Normalized CD8 area was determined by taking the area staining positive for CD8 normalized to the entire area of the tissue section determined by autofluorescence. Models shown in b are indicated. Spearman correlation coefficient (r) and corresponding P-value indicated above plot. Error bars represent sample standard errors. d Correlation of ATM pS1981 determined by RPPA with gene expression values for CCL5, CXCL9, and CXCL10 in the MDST model (N = 11 MDST models). Spearman correlation coefficient (r) and corresponding P-value indicated above plot. Error bars represent sample standard errors

Fig. 8

Neoantigen-independent CTL infiltration correlates with increased ATM across cancer lineages. a Spearman correlation coefficient between CTL score and total ATM determined by RPPA (blue bars), as well as between CTL score and neoantigen levels (red bars). Individual plots shown in Supplementary Fig. 8A. P-values determined from Spearman correlation coefficient and adjusted by Holm-Sidak method. *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 1 × 10^-4. b Correlation coefficient of total ATM protein and CTL levels. Cancers that had significant positive correlation with neoantigen levels showed a lower ATM/CTL correlation than those that did not (Non-neoant. Corr). P-value determined by Wilcoxon rank-sum test. Lines represent median and interquartile range. c Determination of correlation coefficients in M-class and C-class subsets. Correlation coefficients between CTLs and neoantigen load or ATM protein level were assessed for the bulk population (All), as well as M-class and C-class subsets. See Supplementary Fig. 9