B cell-reactive neoantigens boost antitumor immunity

Abstract

B cell involvement in neoantigen-driven antitumor immunity remains largely unexplored because of challenges in predicting B cell responses. Here, we developed a method to identify B cell epitopes by characterizing >437,000 peptides tested for IgG binding and >370 million B cell receptor (BCR) clones. Our single-cell BCR sequencing of pre- and post-severe acute respiratory syndrome coronavirus 2 vaccination validates the performance of this method. Mouse vaccination experiments demonstrate that B cell neoepitopes enhance immune responses, driving BCR expansion and tumor regression. Genomic analysis of >8000 The Cancer Genome Atlas (TCGA) samples reveals an inverse correlation between predicted B cell reactivity and mutation allele frequencies, indicating B cell-mediated neoantigen elimination. Applying our multiomics model to checkpoint blockade responses in 2074 patients highlights the clinical relevance of B cell neoepitope prediction. A meta-analysis of 11 personalized vaccine trials involving 1739 neoantigens suggests that incorporating B cell neoepitopes may improve vaccination efficacy. These results underscore the significance of B cell-reactive neoantigens in antitumor immunity.

Fig. 1.. Construction and evaluation of the B cell epitope prediction model.

(A) Schematic illustration of our model that identifies B cell epitopes, DeepNeo-BCR, on the basis of the CNN architecture. The input matrix contains preferential binding strength between amino acids, and SHM frequency information is used to weight the IGHV residues that are prone to SHMs. The CNN models that predict binding between peptides and an IGHV allele were built independently and later combined into an ensemble model using linear regression. AA, amino acid. (B) Comparison of single-allele models, smaller-scale ensemble models, and larger-scale ensemble models in terms of ROC-AUC, PR-AUC, F1 score, and accuracy score. (C) ROC-AUC curve of validation and test sets. The validation AUC curve was plotted as the mean of threefold cross-validation. (D) PR-AUC curve of validation and test sets. The validation curve was plotted as the mean of threefold cross-validation. (E) Comparison of DeepNeo-BCR with previous linear B cell epitope prediction tools. Shown are the ROC-AUC values reported by each publication.

Fig. 2.. Validation using SARS-CoV-2 vaccination data.

(A) Schematic illustration of identifying BCR clones expanded by vaccination from single-cell BCR sequencing (BCR-seq) data. Pre- and postvaccination BCR sequences of each patient sample were clustered according to sequence similarity. The clones were divided into expanded clones and nonexpanded clones on the basis of clonal size change. (B) Odds ratio (OR) of BCR clones with high DeepNeo-BCR scores having undergone clonal expansion as determined by the difference in cell fraction between pre- and postvaccination. The P value of the odds ratio was calculated with the Fisher exact test. (C) Odds ratio of the high-scoring BCR clones belonging to the group of larger (greater than the mean cell fraction) clones using the postvaccination data. The odds ratio P value was calculated with the Fisher exact test. (D) Sequence logo of the three expanded clones with the highest DeepNeo-BCR scores. (E) Cell type annotation of 65,542 cells across 24 samples after vaccination. Cell types with >100 cells are marked on the UMAP plot. (F) Difference in the mean DeepNeo-BCR score between memory and naïve B cells. The statistics were calculated using the Welch two-sample t test. (G) Odds ratio of the high-scoring BCR clones belonging to the memory B cell population. The odds ratio P value was calculated with the Fisher exact test. (H) Schematic illustration of finding IgG-reactive hotspot epitopes from Li et al. (29). The sequence of BNT162b2, which is part of spike protein, were split into epitopes and were tested for B cell reactivity by adding sera of patients infected with SARS-CoV-2. (I) Distribution of DeepNeo-BCR prediction scores against hotspot and nonhotspot epitopes for expanded (top) and nonexpanded BCR clones (bottom). The D statistic was calculated using Kolmogorov’s D statistic.

Fig. 3.. Effect of B cell reactivity in mouse neoantigen vaccination.

(A) IFN-γ signals were measured for TCR(+)BCR(high) and TCR(+)BCR(low) vaccines (n = 44 peptides and 264 ELISpot assays). The BCR-high and BCR-low vaccines were defined according to the maximal difference between the mutant and wild-type DeepNeo-BCR scores of each potential neoepitope. Individual ELISpot results were compared between the TCR(+)BCR(high) and TCR(+)BCR(low) vaccines. For each vaccine, the immunogenicity assay was conducted in duplicate for each of three mice to ensure robustness. Colored dots indicate individual mice in which the vaccine was tested. (B) Tumor volumes (mm³) were measured longitudinally for the vehicle (n = 17), B16F10-1 (n = 18), and B16F10-2 (n = 18) groups. (C) To directly assess the role of B cells in antitumor immunity, B cell depletion was performed by using anti-CD20. Tumor volumes (mm³) were measured longitudinally for the vehicle (n = 9), B16F10-2 (n = 10), and B16F10-2 + anti-CD20 (n = 10) groups. (B and C) P values were calculated using the two-tailed t test. (D) Comparison of tumor volumes at day 17 between B16F10-2 (n = 10) and B16F10-2 + anti-CD20 (n = 10). Statistical significance was assessed using the one-tailed t test. (E) Correlation plot of Shannon diversity and tumor size on the day of sample collection (day 14) for 34 mice subjected to BCR sequencing (vehicle, n = 11; B16F10-1, n = 13; B16F10-2, n = 10). Statistical analysis was conducted with the Pearson correlation coefficient. (F) IGHV gene usage across vaccine treatment groups. Top 18 IGHV genes with the highest usage across three groups are shown. (G) Morisita-Horn overlap of the BCR repertoire among mice that showed tumor regression when treated with different neoantigen vaccines. Each group is annotated and highlighted with a black box. (H) Alignment of representative CDR3 regions in each treatment group’s regressive tumors.

Fig. 4.. Analysis of B cell reactivity in TCGA pan-cancer data.

(A) Magnitude of neoantigen depletion signal for TCR(+)BCR(+) and TCR(+)BCR(−) epitopes according to immune phenotypes. For each immune phenotype, TCGA samples were divided into high-scoring and low-scoring groups. The magnitude of neoantigen depletion, as measured by the difference in the mean VAF of TCR(+)BCR(+) or TCR(+)BCR(−) mutations between high-immune versus low-immune samples, is plotted for each immune phenotype. (B) Cancer type–wise correlation between BCR clonality and BCE load. (C) Cancer type–wise correlation between BCR entropy and BCE load. (D) Correlation between BCE or TCE load and immune repertoire diversity measurements. All correlation analysis was performed with the Pearson correlation test.

Fig. 5.. Construction of machine learning models for predicting the response to ICB therapy.

(A) Cohort information regarding cancer type, ICB type, sampling time, and clinical response. The heatmap illustrates the availability of omics data categorized by cancer types. The samples without available data were denoted by the gray color. The stacked bar plot presents the proportion of samples with clinical benefit across cancer types. DCB, durable clinical benefit; NCB, no clinical benefit. Created in BioRender. An, J. (2025) https://biorender.com/psv8ito. (B) Model performance evaluated by ROC-AUC in the training (n = 1659) and test (n = 415) sets. The gray dashed line indicates y = x. (C) Box plot illustrating the differences in predictive probabilities across models between DCB and NCB patients in the training (n = 1659) and test (n = 415) sets. Nominal P values were calculated by the two-sided Wilcoxon signed-rank test. The box plot displays the 25th and 75th percentiles with the median represented by the center bar, and the whiskers show the farthest outliers within 1.5 times the interquartile range. (D) Confusion matrix of the XGBoost model in the training (n = 1659) and test (n = 415) sets between true (x axis) and prediction (y axis) results. (E) Kaplan-Meier survival probabilities for overall survival events between the groups predicted by the XGBoost model as DCB and NCB in the training (n = 1659) and test (n = 415) sets. The P values were obtained by the two-sided survival test.

Fig. 6.. Contribution of BCE load to explaining the therapeutic response to ICB.

(A) Bar plot indicating the feature importance measures calculated from each model (XGBoost, HistGBM, and LightGBM). (B) Box plot displaying the distribution of two-sided nominal P values resulting from the bootstrapping test of survival analyses conducted 5000 times, comparing neoepitope load with TMB. Loads of TCE class 1 and class 2 and BCE were detected by DeepNeo-TCR and DeepNeo-BCR, respectively. The gray dashed line indicates P = 0.05. The box plot displays the 25th and 75th percentiles with the median represented by the center bar, and the whiskers show the farthest outliers within 1.5 times the interquartile range. (C) Forest plot showing the odds ratios calculated through logistic regression for response in a meta-analysis conducted across the multiple ICB cohorts. The dashed line signifies an odds ratio of 1. The accompanying bar plot illustrates the nominal two-sided P values obtained from the meta-analysis result of the logistic regression analysis. The dashed line on the bar plot marks the significance threshold of P = 0.05. Features with P values less than 0.05 are highlighted in orange. (D) Forest plot showing the hazard ratios calculated through coxph regression for overall survival in a meta-analysis conducted across the multiple ICB cohorts. The dashed line on the forest plot signifies a hazard ratio of 1. The accompanying bar plot illustrates the nominal two-sided P values obtained from the meta-analysis result of the coxph regression analysis. The dashed line on the bar plot marks the significance threshold of P = 0.05. Features with P values less than 0.05 are highlighted in orange.

Fig. 7.. Effect of B cell reactivity in human personalized neoantigen vaccination.

(A) Concept figure illustrating a vaccine with the somatic mutation marked in red, potential neoepitopes, and predicted BCR reactivity. From a vaccine sequence, multiple 15-nucleotide oligomer peptides were retrieved and tested for B cell reactivity using DeepNeo-BCR. (B) Correlation between the DeepNeo-BCR score and response to vaccines at the patient level. The highest DeepNeo-BCR score among all neoepitopes per patient and patient-level immune response [% of IFN-γ(+) vaccines] were subjected to the Pearson correlation test. (C) Example of a patient case with a high immune response rate (top) and one with a low immune response rate (bottom). Shown are the DeepNeo-BCR scores for potential neoepitopes derived from each vaccine. The scores are shown in the order of neoepitopes that can be obtained by sliding through the vaccine sequence, as illustrated in (A). The red palette indicates IFN-γ(+) vaccines, and the blue palette indicates IFN-γ(−) vaccines. (D) Odds of a patient that received a high fraction of B cell–reactive neoepitopes (>20% BCEs) showing a strong immune response (>70% of delivered vaccines). The odds ratio P value was calculated with the Fisher exact test.