Netie: inferring the evolution of neoantigen-T cell interactions in tumors

Abstract

Neoantigens are the key targets of antitumor immune responses from cytotoxic T cells and play a critical role in affecting tumor progressions and immunotherapy treatment responses. However, little is known about how the interaction between neoantigens and T cells ultimately affects the evolution of cancerous masses. Here, we develop a hierarchical Bayesian model, named neoantigen-T cell interaction estimation (netie) to infer the history of neoantigen-CD8⁺ T cell interactions in tumors. Netie was systematically validated and applied to examine the molecular patterns of 3,219 tumors, compiled from a panel of 18 cancer types. We showed that tumors with an increase in immune selection pressure over time are associated with T cells that have an activation-related expression signature. We also identified a subset of exhausted cytotoxic T cells postimmunotherapy associated with tumor clones that newly arise after treatment. These analyses demonstrate how netie enables the interrogation of the relationship between individual neoantigen repertoires and the tumor molecular profiles. We found that a T cell inflammation gene expression profile (TIGEP) is more predictive of patient outcomes in the tumors with an increase in immune pressure over time, which reveals a curious synergy between T cells and neoantigen distributions. Overall, we provide a new tool that is capable of revealing the imprints left by neoantigens during each tumor's developmental process and of predicting how tumors will progress under further pressure of the host's immune system.

Extended Data Fig. 1

Applying netie on the simulation data. Two more simulation datasets (a) and (b,c). (a) Simulation setting of the second dataset and netie’s inference results. (b) Simulation setting of the third dataset. (c) Netie’s inference results for the third dataset. The same simulation and analysis procedures, as in Fig. 1e–g, were carried out.

Extended Data Fig. 2

Immune selection pressure variations correlate with the genotypes of the tumor and tumor clones. (a) The top genes with smallest Wilcoxon test P values comparing the immune pressure variations in the tumor clones with and without mutations in each gene. (b) Boxplots of the immune pressure variation (ac) in the tumor clones with and without mutations in SETDB1 and FN1. (c) Enriched GO terms of the genes with Wilcoxon test P value<0.05. For boxplot in (b), box boundaries represent interquantile ranges, whiskers extend to the most extreme data point which is no more than 1.5 times the interquartile range, and the line in the middle of the box represents the median.

Extended Data Fig. 3

The immune selection pressure scores of the shared and the private tumor clones of the other MDACC lung cancer patients with multi-region sampling. The center of the error bar represents the inferred immune pressure variation (ac), and the error bars represent 95% confidence intervals.

Extended Data Fig. 4

Further validating the implication of Netie classifications for prognosis of patients (a,b) Multivariate analysis testing the association between TIGEP and overall survival in INisp and DEisp patients. (a) INisp patients, (b) DEisp patients. The association between TIGEP and overall survival is tested in a CoxPH model, with multi-variate adjustment for pathological stage, gender, age, and tumor types.

Extended Data Fig. 5

Analyses as conducted in Fig. 4a, but with the clonality inference conducted by PhyloWGS and SciClone (a), or with the neoantigens predicted by MHCflurry (b). The TCGA LUAD cohort was employed as an example.

Fig. 1

The rationale for inferring the history of anti-neoantigen immune pressure during tumor development. (a) Clonal composition of a hypothetical tumor. The circles refer to the proportions of tumor cells carrying each variant. The circles are colored according to the clones, to which they are assigned to. The histogram shows the distribution of the cellular prevalence of the variants. (b) Phylogenetic relationship between clones. All clones, inferred in (a), are derived either from the normal tissue or from a parent clone. Clone 1 breaks into two clones as two tumor cells are born (Clones 2 and 3). (c) Evolution of mutations and their prevalences within clones. The different shades of pink refer to nesting clones. (d) Inferring immune selection pressure from neoantigens. The blue curly shapes refer to neoantigens associated with each mutation (circle). (e) The setup of the simulation data, where the assumed clones and their parental relationships were shown. (f) The posterior density curves of the random variables to be estimated, with the 95% highest posterior density intervals presented by blue bars on the x-axes. The vertical red lines are located at the true assumed values. (g) Trace plots showing the convergence of the netie estimates of the random variables around the true values, throughout the MCMC iterations. (h) The potential scale reducing factors (PSRFs) for all the inferred variables of the simulation dataset in Fig. 1e. “ac” is the inferred trend of change in anti-tumor selection pressure for each clone. “bc” and “pi” are the posterior estimates of the other variables in the Bayesian model.

Fig. 2

Immune selection pressure variations correlate with the phenotypes of the tumor and tumor clones. (a) Applying netie on the TCGA plus the kidney cancer data. The percentages of the patients with high “a” (a >0 in more than 70% MCMC iterations) and low “a” (a<0 in more than 70% iterations) are shown for each tumor type. (b) Circos plots showing the enriched pathways in the genes that are differentially expressed between INisp and DEisp patients (a> or <0 in more than 50% of iterations). Left: KIRC; right: SKCM. Only the top pathways are shown in each panel for ease of presentation. (c) The number of enriched immune-related pathways found in the genes differentially expressed between INisp and DEisp patients, for each cancer type. (d) The top differentially enriched pathways between INisp and DEisp patients, detected by GSEA. For this analysis, all patients regardless of cancer types were combined. The GSEA test was applied for the calculation of P values. (e) Volcano plot showing the genes that are differentially expressed between INisp and DEisp patients of SARC. A positive value on the X axis means the gene is up-regulated in the INisp patients. Two-sided T-test was applied. (f) A heatmap showing the differential expression of HAVCR2, LAG3, IL-2, IFNG, and TNF, in all cancer types. Orange refers to higher expression in INisp patients, and blue refers to higher expression in DEisp patients.

Fig. 3

Netie is capable of performing multi-sample joint analyses. (a) Netie analysis of the multi-site samples of one MDACC lung cancer patient (Patient ID 886403). The tumor clones were visualized in the phylogenetic tree plot and the fish plot. (b) The immune selection pressure scores of the shared (N=1) and the private tumor clones (N=12) of this patient in (a). (c) Netie analysis of the pre-treatment and post-treatment samples from the Riaz cohort. The immune selection pressure scores were also visualized in barplots for comparison between clones that occurred in pre-treatment samples and new clones that occurred only in the post-treatment samples. The numbers of clones shown in the barplots are as follows. Pt3: N(Pre)=2, N(Post)=1. Pt9: N(Pre)=2, N(Post)=2. Pt10: N(Pre)=3, N(Post)=1. Pt11: N(Pre)=1, N(Post)=1. Pt26: N(Pre)=1, N(Post)=1. Pt27: N(Pre)=1, N(Post)=1. Pt31: N(Pre)=1, N(Post)=1. Pt89: N(Pre)=3, N(Post)=1. (d) Boxplots of the expression levels of the T cell exhaustion signature, comparing the pre-treatment and post-treatment samples. N(pre-treatment)=8 and N(post-treatment)=8. (e) GO analysis of the genes differentially expressed between the pre-treatment and post-treatment samples. The lengths of the bars are proportional to the −log(P value) of the GO analysis. The GSEA test was applied for statistical test. For barplots in (b)-(c), the center of the error bar represents the immune pressure variation (ac) and the error bars represent 95% confidence intervals. For boxplot in (d), box boundaries represent interquantile ranges, whiskers extend to the most extreme data point which is no more than 1.5 times the interquartile range, and the line in the middle of the box represents the median.

Fig. 4

History of immune selection pressure is predictive of tumor prognosis and response to immunotherapy. (a-b) The patients were dichotomized into the INisp and DEisp groups. Survival analyses were performed by examining the association between TIGEP levels (top 50% vs. bottom 50%) and overall survival in each group. (a) all patients; (b) patients from immunogenic cancer types. (c,e) The immunotherapy-treated patients were dichotomized into INisp and DEisp groups, and each group is further split into patients with high and low TIGEPs (50% vs. 50% split). The proportion of responders in each group (INisp/DEisp, high/low T cell) was examined. (c) all patients, (e) patients from immunogenic cancer types. (d,f) The odds ratios of the enrichment of responders and patients with high TIGEPs was calculated and compared between INisp and DEisp patients. The “high TIGEP” patients were selected based on a number of cutoffs to assess the robustness of this analysis. (d) all patients (upper 30% TIGEP level as cutoff: Pval=0.0048, 40%: Pval=0.0040, 50%: Pval=0.0056, 60%: Pvall=0.0099) (f) patients from immunogenic cancer types (30%: Pval=0.032, 40%: Pval=0.027, 50%: Pval=0.016, 60%: Pval=0.047). For barplots in (d) and (f), the center of the error bars represents the odds ratio and the error bars represent 95% confidence intervals. One-tail p values of CoxPH model were shown in (a-b). One-tail P values of Chi-squared test were shown in (c, e).