Molecular and genetic properties of tumors associated with local immune cytolytic activity

Abstract

How the genomic landscape of a tumor shapes and is shaped by anti-tumor immunity has not been systematically explored. Using large-scale genomic data sets of solid tissue tumor biopsies, we quantified the cytolytic activity of the local immune infiltrate and identified associated properties across 18 tumor types. The number of predicted MHC Class I-associated neoantigens was correlated with cytolytic activity and was lower than expected in colorectal and other tumors, suggesting immune-mediated elimination. We identified recurrently mutated genes that showed positive association with cytolytic activity, including beta-2-microglobulin (B2M), HLA-A, -B and -C and Caspase 8 (CASP8), highlighting loss of antigen presentation and blockade of extrinsic apoptosis as key strategies of resistance to cytolytic activity. Genetic amplifications were also associated with high cytolytic activity, including immunosuppressive factors such as PDL1/2 and ALOX12B/15B. Our genetic findings thus provide evidence for immunoediting in tumors and uncover mechanisms of tumor-intrinsic resistance to cytolytic activity.

Figure 1. Immune cytolytic activity (CYT) varies across tumor types and is associated with suppressive factors

(A) Cytolytic activity (CYT), defined as the log-average (geometric mean) of GZMA and PRF1 expression in transcripts per million (TPM), is shown for each of 18 TCGA tumor types and normal tissues. Normal tissue samples include TCGA controls and GTEx samples, excluding smokers for lung tissues. Boxes in box plot represent interquartile ranges and vertical lines represent 5^th–95^th percentile ranges, with a notch for the median. P-values are unadjusted and calculated by Wilcoxon rank-sum test (comparison to relevant normal), and asterisks denote events significant at 10% FDR. (B) The correlation of a gene with CYT across all tumor types is shown (X-axis) relative to its relative expression in CTL/NK cells. Top right, genes expressed in CTL/NK cells that are associated with CYT. Bottom right, non-CTL/NK genes associated with CYT. Average Spearman correlation of expression with CYT was calculated across 18 tumor types. Y-axis: for each gene, median expression in NKs and CTLs divided by median expression in non-hematopoietic cells using CAGE data from Fantom5. See also Data S1 and Table S1.

Figure 2. Viral infection is tumor-specific and associated with higher CYT in a subset of tumor types

(A) Rates of viral infection, as defined by viral RNA-Seq read counts exceeding those observed in GTEx, for tumor types exhibiting at least one case. Isolated cases of several other viruses were also observed. (B) Distribution of CYT in tumor samples with (+) or without (−) viral infection. In tumor types affected by multiple viruses, “negative” samples include only those negative for all viruses. Box plots as in Figure 1. P-values are according to Wilcoxon rank-sum test. See also Data S2 and Table S2.

Figure 3. Count of predicted antigenic mutations per sample is linked with cytolytic activity and selectively depleted in certain tumor types

(A) Local regression curves showing significant relationships between CYT and total mutation count in eight tumor types (p<0.1, Spearman rank correlation), plus melanoma (dotted line). Curves span the 5^th to 95^th percentile of the mutation count variable. Colors correspond to tumor type and are the same as appear in Figure 1. (B) Analogous to (a), but based on the count of point mutations predicted to yield an antigenic neo-epitope. Potential for antigenicity was defined based on gene expression and potential to bind the corresponding patient’s imputed HLA with high affinity. (C) For each tumor, the count of point mutations predicted to generate neo-epitopes was divided by the total count of non-silent point mutations to yield B_obs/N_obs. This observed ratio was compared to an expected ratio, B_pred/N_pred, estimated from the mutational spectra of the silent point mutations in the given sample using an empirical model (Methods). The ratio of the observed and predicted ratios represents the relative deviation of the neo-epitope rate from expectation. P-values reflect Wilcoxon rank-sum tests for deviation from 1. (D) Analogous to Figure 3C, but using neo-epitope prediction based on randomly re-permuted HLA genotype assignments (across patients). Asterisks denote trends significant at 10% FDR for all panels. See also Data S3, Table S3, and Table S4.

Figure 4. Endogenous retroviruses tied to local immunity

(A) RNA-Seq-derived ERV expression in reads per million (RPM) across 18 TCGA tumor types and 27 non-tumor tissue types (from TCGA and GTEX) for three elements found to be tumor-specific. The expression ranges (minimum value to maximum value) are highlighted in orange (for tumor tissues) or green (for non-tumor tissues). (B) Spearman-rank correlations between CYT and ERV expression. Gray squares indicate non-significant association (unadjusted p>0.05) and blank squares indicate no over-expression of the given ERV in the given tumor type (expression strictly below the normal tissue maximum). Asterisks (*) denote Bonferroni-significant associations (adj. p<0.05). See also Data S4 and Table S5.

Figure 5. Gene mutations associated with high or low immune cytolytic activity

(A) Only genes showing pan-cancer significance (adj. p<0.1, red for positive, blue for negative and grey for non-significant association) for non-silent mutation association with CYT are shown in top row. Additional rows, clustered by similarity, show independent significant (unadjusted p<0.05) enrichment upon sub-analysis. The black wedges represent the share of samples exhibiting mutation. Bar plot indicates unadjusted pan-cancer p-values for mutational association with CYT, dashed lines indicating thresholds yielding 1% and 10% FDRs. (B) Association between CASP8 mutational status and FASLG (left axis) and TRAIL (right axis) gene expression (TPM) for tumor types demonstrating at least 5 instances of nonsynonymous CASP8 mutation. Light and dark bars correspond to wild type and (nonsynonymous) mutant samples, respectively. Box plots as in Figure 1. P-values are calculated by Wilcoxon rank-sum test. See also Data S5 and Table S6.

Figure 6. Amplifications and deletions are associated with cytolytic activity in tumors

(A) The significance of association between CYT and amplification (orange) and between CYT and deletion (green) for all genic loci. Rightward lines show unadjusted p values for instances in which the lesion was positively associated with CYT, and leftward lines show unadjusted p values for instances in which the lesion was negatively associated with CYT. Dotted lines represent the significance cutoff yielding 1% and 10% FDRs (and also appear in parts B–E). Labels on the right side mark events significant at the 10% FDR, plus B2M. Potential driver genes appear in parentheses. (B) Locus zoom on the 9p24.2-p23 amplification, each bar corresponding to a single gene. Labeled genes include those with driver potential or those on the locus boundary. (C) Locus zoom on the region containing B2M, which was not genome-wide significant. (D) Locus zoom on the 17p13.1 amplification. (E) Significant associations between CNAs and CYT on the pan-cancer and cancer-specific level (as in Figure 5). Pan-cancer significance was defined at a 10% FDR, and significance for individual tumor types was defined at unadjusted p<0.05. Positive association is indicated with red circles, negative with blue circles, and non-association with gray circles. Black wedges indicate the share of samples exhibiting the event (ie. non-zero GISTIC score at the locus). Bar plot indicates unadjusted pan-cancer p-values for CNAs, sorted by significance, with dashed lines indicating thresholds yielding 1% and 10% FDRs. See also Data S6 and Table S7.

Figure 7. Proposed model for evolution of tumor-immune associations

(A) As the tumor develops, we propose that intrinsic tumor factors – such as mutated neoantigens or viruses – induce local immune infiltrates (blue circles) that include cytolytic effector cells (expressing GZMA/PRF1; red circles) that kill tumors (daggers). These factors are expected to be positively correlated with CYT across tumors. (B,C) Under pressure from cytolytic immune cells, subclones with resistance mutations will grow out over time. (B) One subset of these mutations would enable tumors to evade killing, but does not impact the infiltrate, and are positively correlated with CYT (i.e. higher infiltrate samples are enriched for these mutations). (C) Another subset suppresses the immune infiltrate (i.e. lower infiltrate samples are enriched for these mutations), and is negatively correlated with CYT. Notably, p53 mutations and ALOX amplifications were also significantly negatively associated with CD8A, suggesting a reduction in cell numbers and not just activity. See also Data S7 and Table S8.