NeoPrecis: enhancing immunotherapy response prediction through integration of qualified immunogenicity and clonality-aware neoantigen landscapes

Domain: neuroscience
Published: 2026-01-23
Editorial outcome: minor_revision

Abstract
Despite the transformative impact of cancer immunotherapy, the need for improved patient stratification remains critical due to suboptimal response rates. While neoantigens are central to anti-tumor immunity, current metrics, such as tumor mutation burden (TMB), are limited by their neglect of immunogenicity and tumor heterogeneity. Here we present NeoPrecis, a computational framework designed to improve immunotherapy response prediction by refining neoantigen characterization across MHC-I and MHC-II pathways and by integrating tumor clonality information. NeoPrecis features an interpretable T-cell-recognition model that reveals the critical influence of MHC molecules on TCR recognition beyond mere antigen presentation. Benefit HLA alleles, identified through model-driven contribution analysis, exhibit significant predictive power for patient outcomes in immune checkpoint inhibitor treatment (melanoma: p-value = 0.04; NSCLC: p-value = 0.01). NeoPrecis, via its clonality-aware neoantigen landscape feature, improves immunotherapy response prediction in tumor types with varying prevalence of neoantigens, including heterogeneous NSCLC, which retains more subclonal neoantigens due to lower immunoediting pressure. We thus propose NeoPrecis as a comprehensive evaluative framework for neoantigen assessment by incorporating both immunogenicity and tumor clonality, offering insights into the link between the collective quality of neoantigen landscapes and immunotherapy response.

Introduction
Immunotherapy has revolutionized cancer treatment, with personalized cancer vaccines and immune checkpoint inhibitors (ICI) emerging as two major modalities targeting tumor-specific neoantigens1,2. Cancer vaccines work by boosting the immune system through the expansion of T cell precursors specific for these tumor-specific antigens. Since 2017, numerous cancer vaccine trials have been conducted in cancers such as melanoma3–6, glioblastoma7,8, pancreatic cancer9, and others6,10,11. However, response rates remain suboptimal. For instance, Rojas et al. reported an overall response rate of 50%, yet only 11% of target neoantigens successfully induced a T cell response9. These limitations underscore the need for improved strategies to identify and target immunogenic neoantigens more effectively.

Meanwhile, ICIs have demonstrated significant efficacy in certain cancer types, such as melanoma12,13 and non-small cell lung cancer (NSCLC)14,15, by blocking immune checkpoint pathways and reinvigorating CD8+ T cells recognizing tumor-specific antigens, including neoantigens. Despite their transformative potential, ICIs exhibit variable and often suboptimal response rates16,17. For instance, in melanoma, one of the most responsive cancer types, single-agent ICI achieves a response rate of approximately 40%, while combination therapy increases this to around 60%18,19. However, some patients experience life-threatening immune-related adverse events after receiving ICIs20. Tumor mutation burden (TMB), an FDA-approved biomarker for immunotherapy, acts as a proxy for immunogenic neoantigen abundance21. It is known to correlate with response rates across various cancer types16. High-TMB cancers like melanoma show higher response rates, while low-TMB cancers such as sarcoma have lower efficacy. Nevertheless, within individual cancer types, TMB’s predictive power is limited17. Tumor neoantigen burden (TNB), which relates to the count of MHC-presented neoantigens, offers a more refined metric of the tumor’s immune landscape but fails to account for T-cell receptor (TCR) recognition or the clonal distribution of mutations within tumors22.

Neoantigens are novel peptide antigens derived from tumor-specific mutations, playing a critical role in the immune system’s ability to recognize and target tumors. The immunogenicity of neoantigens is highly associated with their abundance and peptide binding affinity for the major histocompatibility complex (MHC), which displays them on the cell surface23. However, these factors alone are insufficient for a sensitive and efficient selection of neoantigens24. The activation of T cells requires TCR recognition of the peptide-MHC (pMHC) complex. Several indirect metrics have been proposed to evaluate TCR recognition, including the agretopicity ratio (the ratio of MHC-binding affinity between mutated and wild-type peptides)23,25 whereby a mutant peptide could be more easily recognized by the TCR repertoire if the wild-type counterpart was not involved in thymic selection, and foreignness (measuring the sequence similarity to foreign immunogenic antigens)26. Direct approaches include ICERFIRE27, DeepNeo28, and PRIME29,30, which predict immunogenicity based on peptide sequence, and cross-reactivity distance31, which measures the distance between wild-type and mutated peptides in terms of TCR recognition. Despite these efforts, the efficacy of TCR recognition metrics remains limited due to the high diversity of TCR repertoires, shaped by genetic recombination and thymic selection32,33. Additionally, MHC allele variability across individuals leads to distinct peptide repertoires, further constraining the generalizability of predictive metrics, underscoring the need for improved methods.

Tumor clonality further shapes the immune response by influencing the distribution of neoantigens across tumor subclones. Clonal neoantigens, which are derived from mutations shared by all tumor cells, can generate more effective clinical responses than subclonal neoantigens, which are confined to a subset of tumor cells34,35. Highly heterogeneous tumors, marked by a greater proportion of subclonal mutations, may evade immune surveillance due to reduced clonal neoantigen exposure36–38. Thus, CD8+ T cells predominantly react to clonal mutations in early-stage NSCLC34, and NSCLC tumors with a high proportion of subclonal copy-number alterations are at higher risk for recurrence or death than those with a low proportion39. Likewise, high tumor heterogeneity correlates with limited ICI response in both mice and humans in melanoma40 and with poor survival outcomes in a meta-analysis of the TCGA cohort38. These findings suggest that suboptimal clinical outcomes may reflect the inability of the immune cells to cope with tumor heterogeneity and underscore the importance of identifying immunogenic neoantigens, as they are not only important biomarkers of ICI response but also indicators of tumor-immune co-evolution through immune selection.

Historically, neoantigen discovery has predominantly focused on the MHC-I pathway41, which activates CD8+ T cells for direct tumor cell killing. This focus reflects the established role of CD8+ T cells as primary cytotoxic effectors in anti-tumor immunity. However, recent studies have expanded our understanding of CD4+ T cells, highlighting their critical roles in both direct tumor killing and CD8+ T cell priming. Bawden et al. revealed direct cytotoxic effects of CD4+ T cells against melanoma42, while Espinosa-Carrasco et al. emphasized their role in sustaining CD8+ T cell responses through interactions with intra-tumor dendritic cells (DC)43. Further supporting this, Pyke et al.44 and Alspach et al.45 demonstrated the contribution of MHC-II neoantigens to anti-tumor immunity. Moreover, Sears et al.46 found that patients with MHC-II-reliant neoantigen profiles, featuring a higher proportion of MHC-II neoantigens, experience significantly longer-lasting clinical benefits. These findings reinforce the critical role of MHC-II neoantigens and highlight the necessity of identifying both MHC-I and MHC-II neoantigens to fully capture the cooperative interplay between CD8+ and CD4+ T cells in tumor immunity.

To encompass these variables, we present NeoPrecis, a computational framework that predicts neoantigen immunogenicity and constructs a clonality-aware neoantigen landscape model to more accurately evaluate patients’ potential immune responses (Fig. 1). NeoPrecis integrates three critical dimensions—neoantigen abundance, MHC presentation, and TCR recognition—into a unified evaluative framework that incorporates both MHC-I and MHC-II pathways. By accounting for MHC allele variability and tumor clonality, NeoPrecis provides a more precise representation of the tumor’s immune landscape that addresses the limitations of current bulk metrics. This comprehensive approach improves the identification of immunogenic neoantigens and enables better patient stratification for immunotherapy.

Results
Characterizing T-cell recognition in immunogenicity prediction
Measures for assessing neoantigen abundance and binding to MHC molecules are well established. In contrast, predicting T-cell recognition is challenging due to the high variability of TCR repertoires and the still-limited availability of TCR sequencing data. To incorporate a model of T cell recognition into our in silico neoantigen prioritization pipeline, we considered two alternative approaches for assessing TCR recognition without TCR sequencing: peptide-based methods and the cross-reactivity distance (CRD) approach. Cross-reactivity, a well-documented feature of TCR-peptide binding, describes the ability of a TCR to bind multiple peptides with similar sequences47. Peptide-based methods, such as PRIME29,30 and DeepNeo28, primarily model the immunogenicity of target peptides. In contrast, CRD, first introduced by Łuksza et al.31, measures the differences between wild-type (WT) and mutated (MT) peptides, aligning more closely with immune mechanisms, as thymic selection acts to constrain TCR diversity to distinguish self from non-self peptides33. However, existing models overlook the role of MHC-binding motifs, which are critical in shaping TCR-pMHC interactions. To address this limitation, we developed a new immunogenicity model for NeoPrecis, NeoPrecis-Immuno, which incorporates MHC-binding motifs into a refined CRD-based framework, enabling a more nuanced characterization of neoantigens.

NeoPrecis-Immuno embeds MHC-binding core regions (9-mer) using trained amino acid embeddings and MHC-binding motif enrichment, projecting each residue into a latent space specific to MHC alleles. Residue-level embeddings are then aggregated using position-weighted pooling to generate peptide embeddings for WT and MT peptides, highlighting the varying importance of core positions in TCR-pMHC interactions. To capture differences between WT and MT peptides, NeoPrecis-Immuno employs a geometric representation that includes origin, direction, and distance. The distance component predominantly reflects CRD and is further refined with peptide sequence information (origin and direction). Finally, NeoPrecis-Immuno integrates CRD with MHC-binding score as a covariate, allowing the model to further learn immunogenicity beyond what is captured by MHC binding alone, to output a probability that represents neoantigen immunogenicity (Fig. 2a, Methods).

We implemented a two-stage training process using TCR-pMHC binding data from IEDB48 and VDJdb49 (Supplementary Data 1) and T-cell assay data from CEDAR50 (Supplementary Data 2). In the first stage, we constructed a cross-reactive peptide triplet dataset from binding data, where each triplet consists of a seed peptide, a cross-reactive peptide (peptide sharing TCR binding with the seed), and a non-cross-reactive peptide (peptide not sharing TCR binding). We curated 1153 MHC-I and 261 MHC-II cross-reactive pairs and sampled 10 non-cross-reactive peptides for each pair, yielding 11,530 MHC-I and 2610 MHC-II triplets for training. The model was trained to minimize the distance between cross-reactive pairs (seed and cross-reactive peptides) and maximize the distance between non-cross-reactive pairs (seed and non-cross-reactive peptides) (Supplementary Fig. 1a, Methods). In the context of a tumor, mutant peptides more similar to wild-type (cross-reactive with the same TCRs) are less likely to be immunogenic because T cells recognizing these self-similar epitopes are typically eliminated during central tolerance. In contrast, more dissimilar mutants (less cross-reactive) are more likely to be immunogenic, as T cells capable of recognizing these neoantigens escape thymic negative selection and remain available in the peripheral repertoire. Notably, the geometric distance calculated by NeoPrecis-Immuno outperformed the BLOSUM62 and PMBEC (an amino acid similarity matrix for peptide-MHC binding)51 distances in differentiating cross-reactive and non-cross-reactive peptide pairs (Fig. 2b, Supplementary Fig. 2a).

In the second stage, we fine-tuned the model using T-cell assay data with T-cell activation labels, enabling it to learn the immunogenicity of target peptides (Methods). We curated 4176 MHC-I and 87 MHC-II samples predicted to bind MHC (NetMHCpan %rank ≤2 for MHC-I and ≤10 for MHC-II) from CEDAR and split them into training (75%), validation (10%), and testing (15%) sets. The resulting immunogenicity model was subsequently validated on both an internal and an external testing set to benchmark its performance.

We compared NeoPrecis-Immuno to an MHC-binding predictor (NetMHCpan52,53) and immunogenicity models (PRIME29,30, ICERFIRE27, and DeepNeo28) (Methods). The internal testing set, isolated from the CEDAR dataset50 before training, is epitope-based with labels for each pMHC, while the external testing set, derived from an independent cancer cohort (NCI dataset, Parkhurst et al.54), is mutation-based with T-cell assay data for each mutation (Supplementary Data 3). For mutation-based predictions, aggregation across MHC alleles was required; we systemically evaluated aggregation strategies later in the Results (Methods). As there were not many MHC-II data points, internal validation was performed exclusively for MHC-I. To ensure fairness, we identified pMHCs not present in the training datasets of any predictors (n=1147, Supplementary Fig. 2b) and excluded pMHCs with unavailable MHC alleles in any predictor. Among the remaining 438 pMHCs (#positives=228), NeoPrecis-Immuno outperformed NetMHCpan and DeepNeo but performed slightly worse than PRIME and ICERFIRE (Fig. 2c, Supplementary Fig. 2c). We also rebalanced samples through bootstrapping to achieve a 1:3 positive-to-negative ratio, matching the test set balance used in PRIME and ICERFIRE benchmarks and better reflecting realistic immunogenicity rates55, which yielded similar performance rankings (Supplementary Fig. 2d). For external validation, NeoPrecis-Immuno achieved the best performance for both MHC-I (n=1089; #positives=36) and MHC-II (n=1189; #positives=33) predictions. For MHC-I, NeoPrecis-Immuno outperformed all other predictors across AUROC, AUPRC, and PPV metrics. For MHC-II, NeoPrecis-Immuno achieved the highest AUROC while performing comparably to DeepNeo in AUPRC and PPV (Fig. 2d, Supplementary Fig. 2e). Although mutations with no available alleles were excluded, for some mutations, not all alleles were supported (<6 alleles for MHC-I; <10 alleles for MHC-II) by ICERFIRE and DeepNeo (Supplementary Fig. 2f), potentially contributing to their lower performance. Overall, these benchmarks demonstrate that NeoPrecis-Immuno is competitive with other immunogenicity predictors.

To assess model generalizability, we evaluated TCR diversity in the training data and performance on unobserved MHC alleles. For TCR diversity, we analyzed CDR3 coverage in the cross-reactive peptide triplet dataset and found that approximately 50% of cross-reactive peptide pairs are recognized by multiple distinct CDR3 sequences, with more than 50% of amino acid substitutions supported by multiple independent CDR3s (Supplementary Fig. 3a). This multi-CDR3 support indicates that the training data captures generalizable cross-reactivity patterns rather than isolated clone-specific interactions. For MHC allele generalization, we stratified mutations in the NCI dataset by the number of unobserved alleles—those absent from training data. The CEDAR test set contained only 2 unobserved alleles each for MHC-I and MHC-II (Supplementary Fig. 3b), necessitating the use of NCI for this evaluation. In NCI, 28% of mutations had unobserved MHC-I alleles, while nearly all mutations had at least one unobserved MHC-II allele (Supplementary Fig. 3c). Although performance decreased in the unobserved allele group, NeoPrecis-Immuno outperformed other predictors for MHC-I and performed comparably for MHC-II (Supplementary Fig. 3d). We note that alleles unobserved in our training data may have been observed in other predictors’ training sets, potentially contributing to their performance in these cases. These results demonstrate that NeoPrecis-Immuno generalizes effectively across unobserved MHC alleles and that the training data represent diverse TCR recognition patterns, though limited training diversity remains a constraint.

Interpretable MHC-binding motif contributions to T-cell recognition
To evaluate the contribution of each component in NeoPrecis-Immuno, we extracted scores from the model for substitution distance (SubDist), substitution distance with position weighting (SubPosDist), geometric distance (GeoDist; SubPosDist with motif enrichment), and cross-reactivity distance (CRD; GeoDist with sigmoid scaling) (Supplementary Table 1, Methods). BLOSUM62 (BLOSUMDist) and PMBEC (PMBECDist) distances were included as baselines for comparison. These distances quantify the difference between WT and MT peptides. All NeoPrecis-Immuno components outperformed the baseline in both MHC-I and MHC-II contexts (Fig. 3a). For MHC-I, performance steadily improved as additional components were added to the model. However, in MHC-II, differences between components were less pronounced, likely due to the smaller sample size. Sigmoid scaling contributed minimally to CRD, as GeoDist and CRD scores were highly correlated (Supplementary Fig. 4a).

Plotting amino acid embeddings on a 2D projection, NeoPrecis-Immuno residue embeddings grouped amino acids based on their chemical properties, with similar residues clustering together (Fig. 3b). Notably, NeoPrecis-Immuno embeddings differed from BLOSUM62 embeddings; for instance, histidine (H) was positioned closer to tyrosine (Y) and phenylalanine (F), likely due to their shared aromatic rings, while tryptophan (W) was positioned further away, possibly reflecting its larger size (Supplementary Fig. 4c). Thus, NeoPrecis-Immuno embeddings obtained through training on TCR-pMHC interaction data differ from BLOSUM62 embeddings, derived from evolutionary conservation more broadly.

Position factors revealed that the third and fifth residue positions were most influential for MHC-I and MHC-II, respectively (Supplementary Fig. 4b). Interestingly, these are not typical anchor positions for MHC binding, suggesting their significance lies in TCR-pMHC interactions. This hypothesis is supported by the cross-reactive triplet dataset, where amino acid substitutions with significant property changes (BLOSUM62 substitution score ≤0) at these positions were rare (Supplementary Fig. 1b).

The motif enrichment step of NeoPrecis-Immuno played a critical role in refining residue embeddings. Using affine transformations (rotation, reflection, and scaling) for interpretation, we analyzed how motif-enriched residue embeddings align with MHC-binding motifs (Methods). Among these transformations, scaling was the primary contributor to differences between wild-type and mutated amino acids, with scaling factors being specific to each position and allele. This specificity likely reflects the variation in MHC-binding motifs across different positions. Positions dominated by the same amino acid exhibited consistent scaling behavior. For example, positions dominated by proline tended to scale the y-axis, while those dominated by arginine scaled the x-axis (Fig. 3c). These refinements adjusted the relative distances between amino acids in the embedding space, improving the ability of NeoPrecis-Immuno to differentiate wild-type and mutated residues for immunogenic neoantigen prediction.

To quantify the impact of motif enrichment, we analyzed the positive-to-negative distance ratio for different metrics, where positive distance indicates the distance between a wild-type peptide and its immunogenic mutated counterpart, while negative distance indicates the distance between a wild-type peptide and its non-immunogenic mutated counterpart. Neoantigens were grouped by allele and position to minimize bias, and the mean distances for positive and negative pairs were used to compute the ratio. GeoDist demonstrated a larger positive-to-negative distance ratio compared to SubDist (Fig. 3d), revealing that motif enrichment enhanced the model’s ability to capture meaningful distance relationships between amino acids and better distinguish immunogenic from non-immunogenic samples.

We further computed an allele benefit score by averaging scaling factors across positions by MHC molecule (Supplementary Data 4, Methods). A higher benefit score indicates greater overall scaling of amino acid distance, reflecting increased immunogenicity and potential benefits for tumor immunity. To analyze trends, we categorized alleles based on the most dominant amino acid at the anchor position of their binding peptides, which contributes most to the allele benefit score. MHC-I alleles associated with glutamic acid (E), lysine (K), or arginine (R) exhibited the highest benefit scores (Supplementary Fig. 4d), aligning with previously reported benefit HLA supertypes, including HLA-B27, which predominantly bind peptides with arginine at the anchor position, and HLA-B44, which favors peptides with glutamic acid56,57.

Multi-dimensional metrics for neoantigen immunogenicity
Beyond TCR recognition, neoantigen abundance and MHC presentation are two major factors contributing to immunogenicity. We defined eight key features across the three dimensions and evaluated these dimensions with the NCI dataset54 (Methods). For neoantigen abundance, DNA allele frequency (AF) serves as a rough approximation of mutation clonality58, while RNA AF and RNA expression levels quantify the expression of mutated peptides. To capture MHC presentation, patient harmonic-mean best rank (PHBR) measures binding affinity across an individual’s MHC alleles59, and robustness indicates the number of binding MHC alleles. For TCR recognition, in addition to the NeoPrecis-Immuno score, we included agretopicity ratio, proposed by Duan et al.25, and foreignness, proposed by Łuksza et al.26. While PHBR, which aggregates NetMHCpan predictions across MHC alleles using a harmonic mean, has shown robust mutation-centric performance for predicting cell surface presentation59, no standard method exists for aggregating TCR recognition metrics across alleles. To address this, we tested several aggregation strategies (Supplementary Fig. 5a, b, Methods); The masked maximum approach, which prioritizes the peptide with the highest recognition score across alleles that meet a minimum binding threshold necessary for presentation, either outperformed or matched the best binding approach across all three recognition metrics when evaluated on T-cell activity (Fig. 4a). Unlike the best binding approach, which focuses solely on the MHC allele with the highest binding prediction, the masked maximum approach considers all potential binding alleles. These results suggest that TCR recognition is best captured by considering all alleles capable of effective presentation, rather than focusing solely on the strongest binding peptide.

We next assessed the association between mutation-centric features and T-cell activity across all missense mutations (n=6952; CD8+ positive=56; CD4+ positive=66). Consistent with established immune mechanisms, MHC-I metrics were primarily associated with CD8+ T-cell activity, while MHC-II metrics correlated with CD4+ T-cell activity (Fig. 4b). Among the metrics, RNA AF, RNA expression, PHBR, and NeoPrecis-Immuno emerged as the strongest predictors of immunogenicity. DNA AF showed particular importance for MHC-I, while agretopicity ratio and foreignness were more influential for MHC-II predictions.

To refine predictions, we integrated the mutation-centric features into a logistic regression model. After performing feature selection, abundance features (A): DNA AF, RNA AF, and RNA expression; presentation feature (P): PHBR; and recognition feature (R): NeoPrecis-Immuno were included in the final integrated multi-dimensional model, NeoPrecis-Integrated (Supplementary Fig. 6a, Methods). Cross-validation repeated 100 times on the NCI dataset demonstrated that the A+P+R and A+R achieved similar performance (average AUROC: 0.848 vs. 0.849 for MHC-I; 0.806 for MHC-II), both substantially outperforming single-feature models (Fig. 4c, Supplementary Fig. 6b). This occurs because NeoPrecis-Immuno already incorporates MHC binding predictions as a covariate, effectively capturing presentation information. Feature importance analysis confirmed that recognition dominates in A+P+R while presentation contributes minimally, whereas presentation becomes important in A+P models (Fig. 4d, Supplementary Fig. 6c). These results demonstrate that integrating multi-dimensional features improves both immunogenicity prediction compared to single features alone, with NeoPrecis-Immuno serving as a comprehensive metric that captures both T-cell recognition and presentation information.

Advancing immunotherapy response analysis with NeoPrecis
TMB and TNB are designed to quantify the number of potentially immunogenic neoantigens, but often underperform due to their limited consideration of immunogenic factors. To enhance neoantigen identification, we applied NeoPrecis to a meta ICI cohort comprising five melanoma60–64 and three non-small cell lung cancer (NSCLC)35,65,66 cohorts (Supplementary Table 2, Methods). This dataset included 695 samples with whole-exome sequencing data. We filtered for patients with pre-treatment biopsies and no prior ICI treatment to avoid pre-existing immune effects (Fig. 5a). One sample that failed clonal analysis and nine samples lacking RECIST labels were excluded, resulting in 525 patients (277 melanoma and 248 NSCLC). For subsequent analyses, ICI response was defined as complete response (CR) and partial response (PR), while non-response was defined as stable disease (SD) and progressive disease (PD). Overall survival (OS) was assessed for melanoma, while progression-free survival (PFS) was evaluated for NSCLC due to data availability (Supplementary Data 5). To assess cohort compatibility, we evaluated TMB distributions and mutation profiles (Methods). Pairwise Mann-Whitney U tests with FDR correction identified significant TMB differences in only a few cohort pairs (Supplementary Fig. 7a). Principal component analysis of binary mutation profiles revealed no obvious cohort-specific clustering, with silhouette scores of -0.36 for melanoma and -0.30 for NSCLC indicating minimal cohort-driven structure (Supplementary Fig. 7b).

First, we assessed the allele-level contribution to ICI treatment outcomes using the allele benefit score, derived from the NeoPrecis-Immuno model, which quantifies how MHC alleles refine neoantigen immunogenicity. A patient-centric allele benefit score (BenefitScore) was calculated as the geometric mean across all alleles (Methods). We then applied a Cox proportional hazards model, incorporating age and sex as confounders, to estimate its prognostic value. BenefitScore-dual, which integrates both MHC-I and MHC-II, was statistically significant in predicting patient outcomes for both melanoma (p-value = 0.04) and NSCLC (p-value = 0.01), and outperformed BenefitScore-I (MHC-I only) and BenefitScore-II (MHC-II only) (Fig. 5b). We further compared its performance against TMB as a baseline. While BenefitScore alone did not surpass TMB, their combination (CombinedScore) improved the overall predictive power (Methods). These findings support that different MHC alleles modulate neoantigen immunogenicity in distinct ways, ultimately influencing ICI treatment outcomes.

Since allele-level assessment does not directly account for mutations, we next evaluated the mutation-level contribution using NeoPrecis-Immuno burden (NPB), defined as the number of mutations passing an NeoPrecis-Immuno threshold of 0.4 (approximately the median) for both MHC-I and MHC-II (Methods). NPB demonstrated superior performance in stratifying immune responses compared to TMB and TNB, as evidenced by both Cliff’s Delta effect size (Fig. 5c) and a non-parametric statistical test (Supplementary Fig. 8a), suggesting that NeoPrecis-Immuno improves neoantigen immunogenicity predictions by accounting for TCR recognition.

To conduct a more comprehensive mutation-level evaluation, we calculated a tumor-centric immunogenicity score, NeoPrecis-LandscapeSum, by summing NeoPrecis-Immuno predictions across all mutations (Methods). For each mutation, the MHC-dual NeoPrecis-Immuno score was derived by multiplying the MHC-I and MHC-II NeoPrecis-Immuno predictions. This multiplication strategy prioritizes mutations capable of activating both CD4+ and CD8+ T-cell responses. Effective anti-tumor immunity requires coordinated CD4+ T-cell, CD8+ T-cell, and dendritic cell responses, which depend on presenting MHC-I and MHC-II neoantigens within the same antigen-presenting cell43. Given intratumoral heterogeneity, a single mutation activating both pathways is more likely to elicit coordinated immunity than relying on separate mutations that may not co-occur in the same cell. We compared multiplication (I×II) with summation (I+II) and maximum (max(I,II)) aggregation methods. While all methods showed comparable performance at the mutation level (Supplementary Fig. 8b), multiplication was selected based on its superior performance at the subclone level (see below, Supplementary Fig. 10). This integrated approach demonstrated superior performance over TMB and TNB in both melanoma (AUROC=0.66) and NSCLC (AUROC=0.75) (Fig. 5d). Given the better performance of MHC-dual multiplication, we refer to the tumor-centric immunogenicity score as MHC-dual multiplication in subsequent analyses. In survival analyses (Methods), patients in the top half of the tumor-centric score showed significantly better outcomes than those in the bottom half (Fig. 5e). NeoPrecis-LandscapeSum outperformed TMB in melanoma and showed comparable performance in NSCLC (Supplementary Fig. 8c).

Lastly, we evaluated the efficacy of the multi-dimensional integrated immunogenicity prediction model (NeoPrecis-Integrated) in patients with RNA-seq data, identifying 187 samples (137 melanoma and 50 NSCLC). The NeoPrecis-Integrated model, trained on the NCI cohort using five selected features (Methods), demonstrated improved performance over NeoPrecis-Immuno in melanoma but not in NSCLC (Supplementary Fig. 9a). Here, the lower purity of NSCLC tumors (Supplementary Fig. 9b) may have impacted the inference of neoantigen abundance metrics, such as DNA AF and RNA AF, implicating the importance of data quality for robust predictions.

Clonality-aware neoantigen landscapes in tumor heterogeneity
Mutation burden or summation metrics cannot account for the distribution of mutations within tumors. To address this limitation, we sought to characterize the tumor-specific neoantigen landscape by incorporating both immunogenicity and tumor subclonal architecture. First, we performed clonal analysis using PyClone67, grouping mutations into clusters annotated with their prevalence in the tumor and mutation number (Methods). We developed two approaches to integrate tumor clonality: NeoPrecis-LandscapeCCF, which calculates a weighted sum of NeoPrecis-Immuno scores using cancer cell fraction (CCF), and NeoPrecis-LandscapeClone, which aggregates NeoPrecis-Immuno scores by cluster architecture (Methods). For the latter, the NeoPrecis-LandscapeClone scores of mutations within a cluster were summed to determine cluster immunogenicity, which was then averaged across clusters weighted by prevalence to derive tumor-centric immunogenicity (Fig. 6a).

We evaluated the same three MHC-I and MHC-II integration strategies (multiplication, summation, maximum) at the subclone level, where subclone-level MHC-I and MHC-II scores were calculated separately and then integrated (Methods). Multiplication consistently outperformed other methods in both melanoma and NSCLC (Supplementary Fig. 10). This superior performance at the subclone level provides additional biological support for the coordination between CD4+ T cell, CD8+ T cell, and dendritic cell. Since tumor cells within the same subclone share mutations, dendritic cells in that subclone can present multiple neoantigens on both MHC-I and MHC-II, enabling coordinated dual-class T-cell responses against the same tumor subclone.

We benchmarked NeoPrecis against existing approaches that also incorporate mutation clonality within tumors, including the Cauchy-Schwarz index of neoantigens (CSiN)68 and the immunoediting-based optimized tumor neoantigen load (ioTNL)69. CSiN profiles the neoantigen landscape by weighting the neoantigen load of each mutation with its variant allele frequency, while ioTNL aggregates weighted neoantigen loads across subclones potentially subject to immune elimination (Methods). Overall, NeoPrecis outperformed other approaches by more effectively qualifying the immunogenicity of neoantigens within subclones, rather than simply quantifying their number (Supplementary Fig. 11a). To ensure this superior performance was not driven by cohort-specific effects, we assessed models within individual cohorts. NeoPrecis-LandscapeClone outperformed TMB in 4 of 5 melanoma cohorts, while NeoPrecis-LandscapeSum better outperformed TMB in the Ravi NSCLC cohort (Supplementary Fig. 11b). When incorporating cohort as a covariate in cross-validation, the relative performance of all approaches remained unchanged (Supplementary Fig. 11c), confirming the superiority of NeoPrecis across cohorts. In melanoma, our clonality-aware methods generally showed better performance, with NeoPrecis-LandscapeClone achieving the best AUROC of 0.69. However, in NSCLC, mutation burden-based methods yielded better results, with NeoPrecis-LandscapeSum demonstrating the top AUROC of 0.75 (Supplementary Fig. 11a).

To investigate the discrepancy between melanoma and NSCLC, we evaluated tumor heterogeneity for each sample. Tumor heterogeneity was quantified using the Gini index, which measures inequality, across both prevalence and size (mutation number) dimensions to better capture the distribution of mutation clusters (Methods). A high Gini index indicates dominance of a few clusters and greater homogeneity. Conversely, a low Gini index reflects a more heterogeneous tumor. Surprisingly, the distributions of prevalence Gini index (pGini) and size Gini index (sGini) were similar between melanoma and NSCLC (Fig. 6b). Tumors were then classified based on their heterogeneity: the first quadrant (high pGini and high sGini) represented more homogeneous tumors, while the third quadrant (low pGini and low sGini) represented more heterogeneous tumors. In NSCLC, patients with heterogeneous tumors exhibited worse outcomes compared to those with homogeneous tumors, although the difference was not statistically significant (Supplementary Fig. 12a). In contrast, no such trend was observed in melanoma. We further analyzed tumor heterogeneity in relation to smoking status in NSCLC. Never smokers had a higher proportion of heterogeneous tumors compared to smokers (Supplementary Fig. 13a). Additionally, they exhibited lower TMB (Supplementary Fig. 13b) and are known to have poorer responses to immunotherapy70, aligning with the worse prognosis observed in the heterogeneous group.

We evaluated the performance of the clonality-aware neoantigen landscape across the heterogeneity-defined groups (Methods, Fig. 6c, Supplementary Fig. 12b). NeoPrecis-LandscapeSum, which does not account for clonality, performed consistently across both groups and tumor types, indicating no subset bias between the groups. In melanoma, NeoPrecis-LandscapeClone consistently outperformed the other methods, regardless of tumor heterogeneity. In NSCLC, however, NeoPrecis-LandscapeClone showed superior performance in the heterogeneous group but underperformed in the homogeneous group. This pattern aligns with the weaker performance of NeoPrecis-LandscapeSum in never smokers (Supplementary Fig. 13c), who tend to have more heterogeneous tumors. These findings suggest that, as predicted, the clonality-aware neoantigen landscape is particularly effective at capturing immunogenic neoantigens in heterogeneous NSCLC tumors.

The clonality-aware neoantigen landscape primarily reflects the impact of clonal neoantigens; its poor performance in homogeneous NSCLC tumors suggests that clonal neoantigens may play a less significant role in these cases. One possible explanation is that these tumors have undergone stronger immunoediting, leading to the selective elimination of highly immunogenic clonal neoantigens. To investigate this, we examined the proportion of clonal mutations (CCF≥0.85)35 and the ratio of MHC-I binding mutations (PHBR≤2) across different smoking statuses (Fig. 6d, e). Smokers exhibited a higher proportion of clonal mutations but a lower fraction of binding mutations, suggesting immunoediting or the intrinsic nature of smoking-induced mutagenesis yielding non-MHC binding peptides. This provides a potential explanation for the weaker performance of the clonality-aware neoantigen landscape in homogeneous NSCLC tumors, as the remaining clonal neoantigens may have survived due to low immunogenicity or immune evasion mechanisms, thereby weakening their predictive value.

Methods
Preparation of the cross-reactive peptide triplet dataset
TCR-peptide binding data were obtained from IEDB48 and VDJdb49. Binding complexes were included if they contained a complete annotation of the beta complementarity-determining region 3 (CDR3), peptide sequence, and MHC allele. To isolate TCR recognition from the effects of MHC binding, peptides with low MHC-binding affinity were removed using a % rank threshold of 2 for MHC-I and 10 for MHC-II, based on predictions from NetMHCpan-4.152 (MHC-I) and NetMHCIIpan-4.353 (MHC-II). After filtering, the dataset included 69,278 complexes for MHC-I and 3165 complexes for MHC-II. Peptide-MHC complexes were then grouped by MHC allele and CDR3 sequence, ensuring that all peptides within a group bound to the same MHC and TCR molecule. Groups containing only one peptide were removed.

To train the immunogenicity model to distinguish mutated peptides from wild-type peptides, we constructed a dataset of cross-reactive peptide triplets. Each triplet consists of a seed peptide (a reference peptide), a cross-reactive peptide (sharing TCR binding with the seed), and a non-cross-reactive peptide (not sharing TCR binding) (Supplementary Fig. 1a). The model was trained to minimize the distance between seed-cross-reactive pairs (similar peptides) while maximizing the distance between seed-non-cross-reactive pairs (dissimilar peptides). To mimic the effect of single amino acid substitutions, the difference between both cross-reactive and non-cross-reactive pairs was constrained to one Hamming distance. Cross-reactive pairs were sampled from each peptide-MHC complex group to ensure functional similarity between the seed and cross-reactive peptides. All possible pairs were sampled, yielding 1153 MHC-I and 261 MHC-II cross-reactive pairs.

For each cross-reactive pair, non-cross-reactive peptides were randomly selected to be one Hamming distance away from the seed but absent from the cross-reactive set, which comprised peptides within the same peptide-MHC complex group. A BLOSUM62 matrix was used as a prior, ensuring that the substitution between seed and non-cross-reactive peptides resulted in a negative BLOSUM62 score, indicating substantial sequence dissimilarity. To enhance dataset diversity, we applied 10-fold generation, ultimately producing 11,530 MHC-I and 2610 MHC-II triplets (Supplementary Data 1).

Preparation of the CEDAR immunogenicity dataset
Cancer-related immunogenicity data were obtained from CEDAR. Samples were included if they contained annotations for the MHC allele, wild-type peptide, mutated peptide, and immunogenicity label (T-cell assay). To ensure accurate modeling of substitutions, samples where the wild-type and mutated peptides had mismatched lengths were removed. To isolate TCR recognition from the effects of MHC binding, peptides with low binding affinity were filtered out using a % rank threshold of 2 for MHC-I and 10 for MHC-II. After filtering, the dataset comprised 4176 MHC-I samples and 87 MHC-II samples (Supplementary Data 2). Then, the dataset was split into training (75%), validation (10%), and testing (15%) sets. The testing set initially contained 640 samples, from which 438 unique peptides—ensured to be unobserved by any of the predictors—were retained for evaluation. Due to the limited number of MHC-II samples, this restriction was not applied to MHC-II.

Preparation of the NCI cohort
The NCI gastrointestinal cancer cohort, derived from Parkhurst et al.54, was utilized to evaluate the immunogenicity model and construct a multi-dimensional integrated model. The cohort comprises 6952 missense mutations across 75 patients, with both MHC-I and MHC-II alleles available for all patients, along with abundance data (DNA allelic fraction, RNA allelic fraction, and RNA expression level). Notably, RNA expression levels were quartile-normalized for each individual. Each mutated peptide underwent T-cell assay testing, resulting in mutation-centric data with labels for each mutation. This approach differs from peptide-centric data (e.g., CEDAR dataset) where each peptide has a label. In total, 56 mutations exhibited CD8+ activity, and 66 showed CD4+ activity. For peptide-centric features, such as MHC-binding affinity and immunogenicity scores, an aggregation approach was applied to obtain mutation-centric predictions. This method will be described in detail in a subsequent section.

To evaluate immunogenicity prediction, mutations were filtered based on abundance and MHC presentation metrics to exclude factors that could bias TCR recognition comparison. Mutations with zero DNA allelic fraction, zero RNA allelic fraction, or no expression level were removed. Additionally, mutations were filtered using a % rank threshold of 2 for MHC-I and 10 for MHC-II. The final dataset for immunogenicity prediction evaluation consisted of 1089 mutations for MHC-I (36 positives) and 1189 mutations for MHC-II (33 positives). For the construction and evaluation of the multi-dimensional integrated model, all 6952 missense mutations were included to account for all three dimensions of the analysis (Supplementary Data 3).

NeoPrecis immunogenicity model
The NeoPrecis immunogenicity model, NeoPrecis-Immuno, is a custom machine learning framework that integrates key immunogenicity determinants, including residue-level embedding, MHC-binding motif enrichment, peptide embedding, geometric representation, and sequence-aware scaling (Fig. 2a).

Input preparation. pMHC complexes are first predicted using an MHC-binding predictor, and the 9-mer core binding region is identified as the input sequence. For each neoantigen, we obtain paired wild-type (WT) and mutant (MT) peptide sequences represented as one-hot encoded matrices: PWT, PMT ∈ RL×20, where L=9 is the core binding region length and each row is a one-hot vector over 20 amino acids.

MHC binding motif representation. To represent MHC allele in a generalizable manner, we constructed binding motifs from NetMHCpan predictions. We randomly sampled 100,000 peptides from the human peptidome (9-mers for MHC-I, 15-mers for MHC-II) and predicted binding affinity for each available allele. Peptides passing the binding rank threshold (2 for MHC-I, 10 for MHC-II) were retained, and their 9-mer binding cores were used to compute position-specific amino acid frequencies. The resulting binding motif M ∈ RL×20 serves as the MHC allele representation, where each position contains the frequency distribution over 20 amino acids.

Amino acid embedding. One-hot encoded amino acids are first transformed using the BLOSUM62 substitution matrix, then projected into a latent space via a single-layer neural network:

EWT = f_embed(PWT) ∈ RL×da (1)

EMT = f_embed(PMT) ∈ RL×da (2)

Emotif = f_embed(M) ∈ RL×da (3)

where da=2 is the amino acid embedding dimension.

Motif enrichment. To capture the complex interactions between peptide sequences and MHC alleles, we introduce a motif-dependent convolutional layer. The model learns K=8 convolutional kernels Wconv ∈ Rda×da, which are modulated by the MHC motif embedding:

Ck = Emotif * Wconvk ∈ RL×da, k=1,...,K (4)

These motif-enriched kernels are then applied position-wise to peptide embeddings:

HWT = Ck ⊙ EWT ∈ RL (5)

HMT = Ck ⊙ EMT ∈ RL (6)

where ⊙ denotes element-wise multiplication followed by summation across the embedding dimension. The K kernel responses are concatenated and projected to the peptide embedding dimension:

EWT_enriched = Wproj [HWT1,...,HWTK] ∈ RL×dp (7)

EMT_enriched = Wproj [HMT1,...,HMTK] ∈ RL×dp (8)

where Wproj ∈ Rdp×K is a learnable projection matrix and dp=2 is the peptide embedding dimension.

Position-weighted peptide pooling. To aggregate position-enriched embeddings into fixed-length peptide representations, we employ position-weighted pooling. Each of the 9 core binding positions has a learnable weight wi (separate parameters for MHC-I and MHC-II):

vWT = ΣL i=1 wi * EWT_enriched,i ∈ Rdp (9)

vMT = ΣL i=1 wi * EMT_enriched,i ∈ Rdp (10)

where E denotes the embedding at position i. This design allows the model to emphasize structurally important positions in the MHC binding groove.

Geometric representation. To capture the distinction between wild-type and mutant peptides in a way that reflects TCR cross-reactivity, we construct a geometric representation in the peptide embedding space. The geometric transformation consists of three components:

vorigin = vWT ∈ Rdp (11)

vdirection = (vMT - vWT) / ||vMT - vWT|| + ε ∈ Rdp (12)

d = ||vMT - vWT||2 ∈ R (13)

where vorigin represents the wild-type peptide position, vdirection is the normalized direction of mutation in the embedding space, d is the Euclidean distance between wild-type and mutant embeddings, and ε is a small constant to prevent division by zero.

Cross-reactivity distance (CRD). To account for sequence context, we compute a scaling factor using a two-layer neural network:

s = f_scale(vorigin; vdirection) (14)

The CRD is then computed as:

CRD = d * σ(s) (15)

where σ is the sigmoid function. Higher CRD values indicate lower predicted TCR cross-reactivity between wild-type and mutant peptides.

Immunogenicity prediction: The final immunogenicity prediction integrates CRD with MHC binding prediction:

y_immuno = σ(Wout * [CRD; MHC binding] + bout) (16)

This produces an immunogenicity score between 0 and 1.

Model training. The NeoPrecis-Immuno model was trained in two stages. In the first stage, it was trained on the cross-reactive peptide triplet dataset using contrastive learning with triplet loss to optimize CRD:

L_triplet = max(0, CRD(a, p) - CRD(a, n) + margin) (17)

where a represents a seed peptide, p is a cross-reactive peptide, n is a non-cross-reactive peptide, and the margin was set to 1. In the second stage, the model was fine-tuned on the CEDAR immunogenicity dataset using binary cross-entropy loss:

L_BCE = -1/N Σi (yi log ^yi + (1-yi) log (1-^yi)) (18)

where yi is the true immunogenicity label and ^yi is the predicted score. During fine-tuning, all layers were frozen, except for the scaling factor and classifier, allowing the model to focus on adjusting immunogenicity predictions. Training was conducted for 100 epochs with a learning rate of 0.005, using the AdamW optimizer85.

Residue embeddings
To interpret NeoPrecis, we utilized three distinct residue embeddings: BLOSUM62 embedding, NeoPrecis-Immuno residue embedding, and motif-enriched embedding. To generate BLOSUM62 embeddings, we applied principal component analysis (PCA) to the BLOSUM62 matrix, extracting the top two principal components as the embedding for each amino acid. NeoPrecis-Immuno residue embedding is also a two-dimensional representation designed to capture amino acid relationships in the context of TCR recognition. The motif-enriched embedding is a position-wise transformation of the NeoPrecis-Immuno residue embedding, incorporating position-specific MHC-binding motifs. Each position-specific motif is represented as a projection matrix, which is then multiplied with the residue embedding to refine its representation based on MHC-binding context. This transformation ensures that amino acid embeddings reflect not only intrinsic residue properties but also their contextual influence within the pMHC complex.

Immunogenicity predictions from other predictors
To benchmark NeoPrecis-Immuno, we evaluated its performance against three immunogenicity predictors: PRIME30, ICERFIRE27, and DeepNeo28. All three models provide MHC-I predictions, but only DeepNeo supports MHC-II. However, due to differences in MHC allele coverage, not all predictors were available for every allele in our datasets. To ensure fair comparisons in the CEDAR internal validation, we restricted the testing set to peptides supported by all predictors. For the NCI external validation, if a predictor did not support certain alleles in a patient, we removed those alleles and computed the mutation-centric immunogenicity score based on the remaining available alleles. Predictions were aggregated across alleles using the masked maximum approach. If a mutation lacked predictions from a given predictor across all its binding alleles, it was excluded from the testing set to maintain consistency in benchmarking.

For PRIME, we used PRIME 2.0, downloaded from GitHub. ICERFIRE was run as a stand-alone application, excluding expression values from its input. DeepNeo predictions were obtained using the DeepNeo-v2 web server, which only supports the beta chain of MHC-II (DQB, DPB). Consequently, alpha chains (DQA, DPA) were ignored when computing mutation-centric DeepNeo scores. These standardized procedures ensured a fair and unbiased evaluation of NeoPrecis-Immuno against existing immunogenicity predictors.

Component contribution analysis
To assess the contribution of individual components in the NeoPrecis-Immuno, we systematically extracted and analyzed key features, including residue embedding, motif enrichment, position weighting, and CRD scaling factors, and compared them to baseline amino acid distance metrics.

For baseline comparison, we used BLOSUM62 and PMBEC51 substitution matrices. BLOSUM62 is an amino acid substitution matrix derived from aligned protein sequences, while PMBEC is specifically derived from peptide-MHC binding data. Since substitution matrices provide discrete scores for amino acid pairs, we applied PCA to each matrix to extract the first two principal components as two-dimensional amino acid embeddings, enabling continuous distance calculations. We then computed distances as the Euclidean distance between wild-type and mutated amino acids in this reduced space (BLOSUM62Dist, PMBECDist).

NeoPrecis-Immuno progressively refines distance calculations through four stages. First, substitution distance (SubDist) uses trained residue embeddings to measure differences between wild-type and mutated peptides. This is then refined into position-weighted substitution distance (SubPosDist), which adjusts SubDist by incorporating positional importance, ensuring that residues within the peptide-binding core are weighted according to their functional contribution to TCR recognition. Further refinement is achieved through geometric distance (GeoDist), which extends SubPosDist by integrating motif enrichment. By incorporating MHC-binding motifs, GeoDist adjusts residue embeddings to better reflect TCR-pMHC interactions. Finally, cross-reactivity distance (CRD) applies a sigmoid scaling function to GeoDist, modifying geometric distances to account for sequence-specific variations in immunogenicity. This scaling ensures that even when geometric distances between wild-type and mutated peptides are similar, differences in sequence composition are considered in immunogenicity predictions.

Affine transformation for the motif-enriched embedding
To interpret how motif enrichment refines residue embeddings, we applied an affine transformation between the original residue embedding and the motif-enriched embedding:

Em = AEr + t (19)

where Em is the motif-enriched embedding, Er is the original residue embedding, A is the transformation matrix, and t is the translation vector. To analyze how this transformation scales the embedding space, we performed singular value decomposition (SVD) on A:

A = UΣVT (20)

where Σ contains the scaling factors, and VT provides the principal directions of transformation. The scaling effect along the original x- and y-axes was computed by projecting Σ onto the original coordinate system using VT:

Sx = |V1x|*σ1 + |V2x|*σ2, Sy = |V1y|*σ1 + |V2y|*σ2 (21)

where Sx and Sy represent the scaling along the x- and y-axes, and v1, v2 are the basis vectors from VT. This analysis quantifies how motif enrichment reshapes residue embeddings, refining their representation in the context of MHC-binding motifs.

Allele benefit scores
Allele benefit scores are designed to evaluate how MHC-binding motifs refine neoantigen immunogenicity by influencing T-cell recognition.

For each allele and position, we calculated scaling factors (Sx and Sy), and the allele benefit score was defined as the average geometric mean of Sx and Sy across all 9-mer positions:

BenefitScore = 1/9 Σi=1 9 sqrt(Sx,i * Sy,i) (22)

A higher allele benefit score results from larger Sx or Sy values, indicating a greater distance between WT and MT peptides, which in turn suggests increased immunogenicity. Additionally, the allele benefit score can be aggregated across individual alleles to provide an overall assessment of an individual’s immune response. The geometric mean was calculated separately for MHC-I and MHC-II alleles, where N is the number of alleles:

MHCBenefitScore = (Πj=1 N BenefitScorej) 1/N (23)

The MHC-dual benefit score was computed as the geometric mean of the MHC-I and MHC-II benefit scores:

MHC-dualBenefitScore = sqrt(MHC-IBenefitScore * MHC-IIBenefitScore) (24)

To integrate with TMB, the combined score was calculated as the product of logTMB and MHC-dual benefit score:

CombinedScore = (logTMB) * (MHC-dualBenefitScore) (25)

Calculation of multi-dimensional metrics
Eight metrics were computed to evaluate neoantigen immunogenicity, categorized into abundance, presentation, and recognition metrics. Abundance metrics included DNA allele frequency (AF), RNA AF, and RNA expression quartile, while presentation metrics comprised robustness and PHBR. Recognition metrics included agretopicity ratio, foreignness, and NeoPrecis-Immuno. Abundance metrics were available in the NCI cohort, whereas for the Ott and ICI cohorts, DNA AF was derived from variant calling results, and RNA AF was calculated using bam-readcount86. RNA expression quartile was computed as the quartile transformation of TPM values obtained from RSEM87 for each individual. Presentation metrics were determined as follows: robustness was defined as the number of binding alleles, using a % rank threshold of 2 for MHC-I and 10 for MHC-II. PHBR (patient harmonic-mean best rank) was computed as the harmonic mean of the % rank of the best-binding peptide for each allele. Recognition metrics included agretopicity ratio, which was calculated as the ratio of MHC binding scores between mutated and wild-type peptides, and foreignness, which measured sequence similarity to foreign immunogenic antigens. Foreignness was computed using the NeoantigenEditing repository (GitHub) by Łuksza et al.31.

T-cell recognition metric aggregation across MHC alleles
To predict mutation-level immunogenicity from allele-specific scores, we evaluated multiple aggregation strategies that account for individual MHC genotypes. Each strategy involves two key decisions: whether to apply allele masking and which aggregation method to use.

Allele masking: Masking excludes alleles predicted not to bind the target peptide, based on NetMHCpan % rank thresholds (≤2 for MHC-I, ≤10 for MHC-II). Let Aall denote all MHC alleles for an individual and Abind ⊆ Aall denote binding alleles. Without masking, all alleles contribute to aggregation; with masking, only Abind is considered.

Aggregation methods: For allele-specific scores sa, we evaluated five aggregation methods.

– Mean: smut = 1/|A| Σa∈A sa

– Harmonic mean: smut = |A| / Σa∈A (sa+ε)-1

– Weighted average: smut = Σa∈A wa * sa, where wa is binding affinity

– Best binding: smut = sa*, where a* = argmin rank(a)

– Maximum: smut = max a∈A sa

We compared all ten combinations (with/without masking × five aggregation methods) on the NCI dataset and found that masked maximum approach performed most consistently across both MHC classes and different immunogenicity metrics (Supplementary Fig. 5b).

To assess the impact of MHC allele aggregation on immunogenicity prediction, we tested various approaches. The first consideration was masking, which excludes alleles that do not bind to the neoantigen based on the % rank threshold of 2 for MHC-I and 10 for MHC-II. The second consideration was the choice of aggregation method, where we evaluated mean, harmonic mean, weighted average (weighted by binding affinity), best binding, and maximum. For example, the masked maximum method selects the highest score among all binding alleles that meet the masking criteria. This approach ensures that only strongly binding alleles contribute to the aggregated immunogenicity score, refining the prediction of neoantigen immunogenicity.

Feature selection for the integrated multi-dimensional model
To identify the most predictive features among the eight immunogenicity metrics, we applied recursive feature elimination with cross-validation (RFECV) using 4-fold cross-validation. Features were eliminated one at a time, and the subset yielding the best cross-validation performance was selected as the final featureset. For consistency, we retained five features for both MHC-I and MHC-II, even though the optimal feature set for MHC-II consisted of only two features. This decision was made because there was no significant difference in performance between using two and five features in MHC-II (Supplementary Fig. 6a), and the smaller sample size for MHC-II limited the reliability of feature selection.

Integrated multi-dimensional model
The integrated model, NeoPrecis-Integrated, is a logistic regression trained on the selected features, with all features undergoing standard normalization based on the training set. The mean and standard deviation calculated from the training set were recorded and subsequently applied to the testing set for consistent normalization. To assess the impact of multi-dimensional integration, we conducted 4-fold cross-validation on the NCI dataset. For ICI response prediction, the integrated model was trained on the entire NCI dataset to ensure robust performance.

Preparation of the ICI cohorts
We compiled eight patient cohorts treated with ICIs, including five melanoma cohorts and three NSCLC cohorts (Supplementary Table 2). In total, 695 patients (443 melanoma and 252 NSCLC) with whole-exome sequencing (WES) data were collected from dbGaP. Raw sequencing data were processed using a standardized pipeline. Variant calling was performed using Nextflow-Sarek (v3.4.2)88,89, with the hg38 reference genome, employing Mutect2 for somatic variant calling, VEP for mutation annotation, and ASCAT for copy number analysis, tumor purity, and ploidy estimation. RNA sequencing data were aligned to hg38 using STAR (v2.7.3a)90 with default settings, and transcript quantification was conducted using RSEM (v1.3.1)87. HLA typing was performed with HLA-HD (v1.7.0)91, ensuring accurate identification of patient-specific MHC alleles. Clonality analysis was conducted using PyClone (v0.13.1)67, with 5,000 iterations and a maximum of 100 clusters, allowing for a detailed assessment of tumor heterogeneity.

Calculation of mutation burden
We computed three distinct mutation burden metrics for comparison: tumor mutation burden (TMB), tumor neoantigen burden (TNB), and NeoPrecis-Immuno burden (NPB). TMB was defined as the total number of non-synonymous mutations in the tumor. TNB accounted for non-synonymous mutations that met an MHC-I binding threshold (PHBR≤2), incorporating MHC-binding constraints into the burden calculation. NPB incorporated both MHC-I and MHC-II binding along with TCR recognition. A mutation was counted if the product of NeoPrecis-Immuno-I (MHC-I) and NeoPrecis-Immuno-II (MHC-II) scores was ≥ 0.16. This threshold was derived from the median NeoPrecis-Immuno score (around 0.4) for both MHC-I and MHC-II, where 0.16 represents the squared value (0.4 × 0.4), ensuring that mutations with strong immunogenic potential were captured.

Calculation of tumor-centric immunogenicity score
Aggregating neoantigen immunogenicity scores into a tumor-centric immunogenicity score is essential for assessing individual immune response. To achieve this, we first computed mutation-centric immunogenicity scores using three distinct metrics: NeoPrecis-Immuno-I, representing the NeoPrecis-Immuno score for MHC-I; NeoPrecis-Immuno-II, representing the NeoPrecis-Immuno score for MHC-II; and NeoPrecis-Immuno-dual, calculated as the product of NeoPrecis-Immuno-I and NeoPrecis-Immuno-II, thereby integrating both MHC-I and MHC-II pathways. Unless otherwise specified, NeoPrecis-Immuno-dual was used as the default metric in our experiments of ICI response prediction.

We explored three methods to aggregate these mutation-centric scores into a comprehensive tumor-centric immunogenicity score. The NeoPrecis-LandscapeSum method involved summing the NeoPrecis-Immuno scores across all identified neoantigens, providing a straightforward cumulative measure. The NeoPrecis-LandscapeCCF approach leveraged the cancer cell fraction (CCF) by applying a weighted sum based on mutation prevalence as inferred from PyClone analysis, thus incorporating the tumor’s cellular composition into the immunogenicity assessment. The NeoPrecis-LandscapeClone method integrated tumor clonality to offer a more nuanced perspective. Using PyClone, mutations were grouped into clusters according to their prevalence in the tumor and mutation size. Within each cluster, NeoPrecis-Immuno scores were summed separately to derive cluster-specific immunogenicity for MHC-I and MHC-II, which were then integrated using multiplication (summation and maximum were also evaluated for comparison). The MHC-integrated, cluster-centric immunogenicity scores were then averaged across clusters, weighted by their prevalence, to calculate the tumor-centric immunogenicity score (Fig. 6a). The three NeoPrecis-Landscape scores were log-transformed to facilitate interpretation within a more manageable range of values.

Clonality-aware approaches for benchmarking
We benchmarked NeoPrecis-Landscape against two other approaches that derive information from tumor clonality, the Cauchy-Schwarz index of neoantigens (CSiN)68 and the immunoediting-based optimized tumor neoantigen load (ioTNL)69.

The CSiN score is calculated by aggregating the normalized product of variant allele frequency (VAF) and neoantigen load across various quality cutoffs for predicted peptide-MHC (pMHC) binding. For each mutation, we considered all possible 8 to 11-mer (MHC-I) and 15-mer (MHC-II) mutant peptides paired with their respective MHC alleles. The neoantigen load (Li) for a given mutation i represents the number of these pMHC pairs that surpass a specific MHC-binding prediction rank percentile cutoff (c). Following the original CSiN methodology, we employed rank percentile cutoffs of 0.375, 0.5, 0.625, 0.75, 1.25, 1.75, and 2. Mutations with a VAF (Vi) below 0.05 were excluded. For samples with over 500 mutations, only the 500 mutations with the highest VAF were retained for analysis. The CSiN score is then computed using the following equation:

CSiN = Σi=1 k log Σn i=1 I(q(i) < c) * Vi * Li(c) (26)

where k is the total number of cutoffs used, n is the number of mutations considered for the patient, and q(i) is the best rank percentile among all pMHCs for mutation i.

The ioTNL score is calculated by summing the product of cancer cell fraction (CCF) and neoantigen load for subclones identified through PyClone analysis67 that exhibit evidence of immune elimination. An immune elimination score (Lj/Sj) is calculated for each subclone to quantify its potential for immune elimination; only subclones with a score below a defined cutoff (e) are included. The ioTNL score is computed using the following equation:

ioTNL = Σj=1 m I(Lj/Sj < e) * Lj * CCFj (27)

where m is the total number of subclones, Lj is the neoantigen load in subclone j, Sj is the number of nonsynonymous mutations in subclone j, and e is the predefined cutoff for the immune elimination score. The selection of the cutoff (e=1.4) was empirically determined to maximize the predictive performance of ioTNL across our datasets. Although this data-driven approach carries a risk of overfitting and potential information leakage, no established biological rationale was available to guide this selection.

Survival analysis
For survival analysis, patients were stratified into two groups based on their target metric scores, with one group comprising the top half and the other comprising the bottom half of the distribution. Due to differences in data availability, we used overall survival (OS) for melanoma and progression-free survival (PFS) for NSCLC. Kaplan-Meier survival curves were generated, and log-rank tests were performed to assess statistical differences between the two groups.

Tumor heterogeneity assessment
To evaluate tumor heterogeneity, we calculated the Gini index at two levels: cluster mutation size and cluster prevalence, capturing diversity in both dimensions. The Gini index quantifies inequality, with a value of 0 indicating perfect equality (all clusters have equal mutation counts or prevalence) and a value of 1 indicating maximum inequality (a single cluster contains all mutations or has a prevalence of 1). The Gini index (G) is computed as follows, where n represents the number of clusters and x denotes either mutation size or prevalence:

G = Σi=1 n Σj=1 n |xi - xj| / 2n * Σi=1 n xi (28)

To classify tumors into homogeneous and heterogeneous groups, we applied a two-dimensional approach based on prevalence Gini index (pGini) and size Gini index (sGini). For each cancer type, tumors with pGini ≤ median and sGini ≤ median were categorized as heterogeneous, while tumors with pGini > median and sGini > median were classified as homogeneous.

Batch effect evaluation
To assess potential batch effects across the combined immunotherapy cohorts, we performed two complementary analyses: distribution-based comparisons and mutation profile analysis.

TMB distribution analysis. For each cancer type, we categorized patients by response status (responder and non-responder). Within each cancer type and response category, we performed pairwise Mann-Whitney U tests using Benjamini-Hochberg adjusted p-values to compare TMB distributions across cohorts. Cohort pairs with FDR-adjusted p-values < 0.05 were considered to have significantly different TMB distributions.

Mutation profile analysis. We constructed a binary mutation matrix (presence/absence) for all protein-coding genes across all patients. Principal component analysis (PCA) was performed separately for melanoma and NSCLC. To quantify cohort-specific clustering, we calculated silhouette scores for the top 10 principal components using cohort labels as the clustering variable. Silhouette scores range from -1 to 1, where values near 0 indicate no meaningful clustering structure, and negative values suggest that cohort assignment does not reflect natural groupings in the data, both supporting the absence of cohort-driven clustering.

Performance evaluation across cohorts. To assess whether batch effects impacted model performance, we conducted within-cohort validation by calculating AUROC for each score separately for each cohort. We then performed 5-fold cross-validation incorporating cohort as a categorical covariate in a logistic regression framework.

Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Data availability
The processed data generated in this study are provided in the Supplementary Data files, including cross-reactive peptide triplets, CEDAR immunogenicity data, processed NCI gastrointestinal cancer cohort data, and the ICI cohort metadata with analysis results. Publicly available datasets analyzed in this study were obtained from the Immune Epitope Database (IEDB, https://www.iedb.org), VDJdb (https://vdjdb.cdr3.net), and the Cancer Epitope Database and Analysis Resource (CEDAR, https://cedar.iedb.org). The mutation-centric immunogenicity data (NCI dataset) were obtained from the supplementary materials of Parkhurst et al. (https://doi.org/10.1158/2159-8290.CD-18-1494). Clinical and genomic data from the reanalyzed ICI cohorts were obtained with the following accession numbers: Hugo et al. (SRP090294 [https://trace.ncbi.nlm.nih.gov/Traces/study/?acc=SRP090294], SRP067938), Van Allen et al. (SRP011540 [https://trace.ncbi.nlm.nih.gov/Traces/study/?acc=SRP011540]), Snyder et al. (SRP072934 [https://trace.ncbi.nlm.nih.gov/Traces/study/?acc=SRP072934]), Riaz et al. (SRP094781 [https://trace.ncbi.nlm.nih.gov/Traces/study/?acc=SRP094781]), Liu et al. (SRP011540 [https://trace.ncbi.nlm.nih.gov/Traces/study/?acc=SRP011540]), Ravi et al. (SRP413932 [https://trace.ncbi.nlm.nih.gov/Traces/study/?acc=SRP413932]), Anagnostou et al. (SRP238904 [https://trace.ncbi.nlm.nih.gov/Traces/study/?acc=SRP238904]), Rizvi et al. (SRP064805 [https://trace.ncbi.nlm.nih.gov/Traces/study/?acc=SRP064805]).

Code availability
Codes for conducting immunogenicity prediction and neoantigen landscape evaluation are deposited at both GitHub (https://github.com/cartercompbio/NeoPrecis) and Zenodo (https://doi.org/10.5281/zenodo.17959604).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Many papers have tried to use HLA-presentation predictors as a proxy for neoantigen identifications for the field of\nimmunotherapy. However, few has explicitly trained for missense mutation recognitions by the immune systems. This paper\nfilled the gap by proposing a machine-learning method and a tailored-training dataset in identifying neoantigens arising from\nmissense-mutations. This also provided a more nuanced method in patient stratification in immunotherapy. My\nrecommendation at this stage is revision.\nMajor comments\n1. Multiple immunotherapy cohorts (5 melanoma cohorts and 3 NSCLC cohorts) were combined to do analyses. How did\nyou adjust for cohort-to-cohort differences? Or sequence batch effects?\n2. The training set that was created was created in mind of finding differences in T cells’ ability to differentiate peptides with\nonly 1 amino acid difference. How many unique T cells were included per pMHC triplet? Are these T cells public clones or\nprivate clones of patients? My concern about this method is whether this can be generalized or is it only accounting for\nspecific T cells that are private a few people?\n3. Training data size for HLA class II is small. How many alleles were you able to train for? How did you account for overfit?\nDo you see any skew towards specific alleles for either HLA class I or II?\nMinor comments:\n1. Please rearrange locations of figures 3B, 3C, 3D to make them read from left to right and top to bottom for easier\nreadability", "reviewer_2": "The authors describe a computational framework that predicts the immunogenicity of neoantigens, with the option to include\ntumor clonality as an additional feature to enhance prediction of response to ICI therapy. The authors show on-par\nperformance of the immunogenicity model when compared with previously published studies using similar approaches.\nHowever, the improvement in performance is not strikingly different and could be highly dependent on the test dataset. I\nbelieve the novelty of the tool could be more focused on the incorporation of clonality information, and how this additional\ninformation can be leveraged to improve predictions across more tumor cohorts. Please see detailed comments below:\nMajor comments:\n1.The authors highlight the use of allele-specific MHC binding motifs as well as position-weighted pooling as novel features\nfor their immunogenicity prediction model. How were the MHC binding motifs and positional weights derived or calculated?\nGiven that motifs are generally derived from peptides that bind or are presented by the respective allele, how does the model\nhandle alleles with no training data and thus no clear binding motif?\n2.Distance based on BLOSUM62 was used throughout the manuscript for multiple comparisons (e.g., comparing to the\ngeometric distance calculated by NeoPrecis-Immuno). Could the authors show similar comparisons using PMBEC, or further\nexplain why PMBEC was not considered a more biologically relevant distance matrix?\n3.For NeoPrecis-Immuno’s evaluation on the CEDAR dataset (internal test set), the positive-to-negative ratio is roughly 1:1\n(line 165). This ratio seems quite different from other benchmark datasets (e.g., test sets used in ICERFIRE and PRIME).\nDoes this ratio reflect the true proportion of peptides that can be recognized by T cells? Additionally, for the external test set,\nthe authors mention aggregation methods across multiple peptides for determining mutation site immunogenicity. Have the\nauthors tried and compared multiple aggregation strategies, and are performance results dependent on the aggregation\nmethod? Do the authors mask positive peptides found in training prior to aggregation and evaluation? Could the authors\nprovide and comment on model performance comparisons at an allele+length-specific level?\n4.Given the relatively small dataset sizes, could the authors clarify the number of alleles and peptide lengths present in the\ntraining/testing data? How reliable is the model when generalizing to alleles not seen during training (do such cases exist in\nthe test set)?\n5.When calculating a tumor-centric immunogenicity score using NeoPrecis-Immuno predictions (line 298), the authors\ndirectly multiply MHC-I and MHC-II predictions for the same mutation. This seems to prioritize a single mutation that\nactivates both a class I and class II response, rather than two separate mutations where one activates class I and the other\nactivates class II. Can the authors elaborate on why they believe this is the appropriate strategy?\n6.The authors primarily focused on demonstrating the incorporation of tumor clonality information for generating a tumor-\nlevel overall score. However, it would be interesting to see whether incorporating clonality at the mutation level (without\nusing DNA VAF as a proxy) could help improve neoantigen immunogenicity predictions.\nMinor comments:\n- When describing the construction of the triplet dataset, it was a bit confusing that the authors describe the positives as non-\nimmunogenic mutated peptides and the negatives as immunogenic (lines 145–147), while in Figure S1A, the seed and\npositive peptides are described as having corresponding TCRs. Given that negative peptides were randomly selected, it is\nunclear why they are described as immunogenic mutated peptides.\n- It would be appreciated if the authors labeled the training dataset sizes for the models (e.g., how many seed peptides were\nused from IEDB/VDJdb, and how many data points in CEDAR were retained for training, given that it was also used for\ninternal validation).\n- For evaluation metrics, could the authors also report PPVn in the context of neoantigen immunogenicity prioritization?\n- How are the core positions for Class II peptides determined?\n- In Figure 3D, could the authors comment on the outlier and why motif-embedding did not improve performance in this\ncase?\n- In Figure 4C/S5B, there doesn’t appear to be a significant gain when comparing A+R with A+P+R. Can the authors\ncomment on why this might be the case? Also, the AUROC values for A+R and A+P+R are extremely close — it would be\nhelpful to have the values labeled directly.\n- In Figure 4D, could the authors also show the feature importance of binding/presentation predictions from NetMHCpan?"}, "author_response": "We thank the reviewer for the helpful comments that have strengthened the manuscript. We have\nconducted more analyses to address the concern about model generalizability and robustness, including\ncohort-to-cohort batch correction in ICI response prediction, CDR3 coverage in the training set, and\nmodel generalizability across unobserved MHC alleles.\nMajor comments:\n1. Multiple immunotherapy cohorts (5 melanoma cohorts and 3 NSCLC cohorts) were combined to do\nanalyses. How did you adjust for cohort-to-cohort differences? Or sequence batch effects?\nThank you for raising this important concern. Although we applied a uniform bioinformatics protocol\nfor somatic mutation calling and HLA typing from raw sequences, we acknowledge that sequencing-\nrelated batch effects between cohorts could potentially influence our findings. To address this, we\nsystematically evaluated potential batch effects across cohorts and assessed their impact on\nimmunotherapy response prediction performance.\nBatch effect:\nSince our study focuses on somatic mutations and immunotherapy response, we evaluated batch\neffects at two levels: tumor mutation burden (TMB) and mutation profiles. For TMB, we performed\npairwise Mann-Whitney U tests with FDR correction to compare TMB distributions across cohorts\nwithin the same cancer type and response category. While most cohort pairs showed no significant\ndifferences, we identified several exceptions, including Riaz vs. Liu and Riaz vs. Van Allen in\nmelanoma non-responders, as well as Ravi vs. Rizvi in NSCLC non-responders (Fig. R1A). For\nmutation profiles, we constructed a binary presence/absence mutation matrix for all coding genes\nacross patients and performed principal component analysis (PCA) separately for each cancer type.\nThough there are some outliers, the top two principal components revealed no obvious cohort-specific\nclustering, with low explained variance (Fig. R1B). Furthermore, silhouette scores for the top 10 PCs\nwere -0.36 for melanoma and -0.30 for NSCLC, indicating minimal cohort-driven structure in the\nmutation data.\nConclusion on Batch effect. Together these analyses suggest that batch effects related to\nnonsynonymous somatic mutations are modest across cohorts.\nBatch correction:\nTo further evaluate whether cohort-related batch effects impact model performance, we conducted\ntwo complementary analyses. First, we assessed performance within individual cohorts and found\nthat NP-LandscapeClone outperformed TMB in 4 out of 5 melanoma cohorts, while NP-\nLandscapeSum outperformed TMB in the Ravi NSCLC cohort (Fig. R1C). Second, we incorporated\ncohort as a covariate in a 5-fold cross-validation framework to directly adjust for potential batch effects.\nThe results remained consistent with our original findings (Fig. S11a), with NP-LandscapeClone\nachieving the best performance in melanoma and NP-LandscapeSum in NSCLC (Fig. R1D).\nIn summary, our batch effect analysis reveals only mild cohort-to-cohort variation, and when\naccounting for this variation, the result remained unchanged with NeoPrecis consistently\noutperforming other methods. We have added these batch effect analyses to the supplementary\nmaterials to provide greater transparency regarding data heterogeneity across cohorts.\nManuscript updates:\nUpdated Figures: Fig. R1A-B has been added as Fig. S7a-b.; Fig. R1C-D has been added as Fig.\nS11b-c.\nUpdated Results:\nLine 315-320\nTo assess cohort compatibility, we evaluated TMB distributions and mutation profiles (Methods). Pairwise\nMann-Whitney U tests with FDR correction identified significant TMB differences in only a few cohort pairs\n(Fig. S7a). Principal component analysis of binary mutation profiles revealed no obvious cohort-specific\nclustering, with silhouette scores of -0.36 for melanoma and -0.30 for NSCLC indicating minimal cohort-\ndriven structure (Fig. S7b).\nLine 393-398\nTo ensure this superior performance was not driven by cohort-specific effects, we assessed models within\nindividual cohorts. NeoPrecis-LandscapeClone outperformed TMB in 4 of 5 melanoma cohorts, while\nNeoPrecis-LandscapeSum outperformed TMB in the Ravi NSCLC cohort (Fig. S11b). When incorporating\ncohort as a covariate in cross-validation, the relative performance of all approaches remained unchanged\n(Fig. S11c), confirming the superiority of NeoPrecis across cohorts.\nUpdated Methods: Line 883-901\nTo assess potential batch effects across the combined immunotherapy cohorts, we performed two\ncomplementary analyses: distribution-based comparisons and mutation profile analysis.\nTMB distribution analysis: For each cancer type, we categorized patients by response status (responder\nand non-responder). Within each cancer type and response category, we performed pairwise Mann-\nWhitney U tests comparing TMB distributions across cohorts. P-values were adjusted using the Benjamini-\nHochberg false discovery rate (FDR) correction. Cohort pairs with FDR-adjusted p-values < 0.05 were\nconsidered to have significantly different TMB distributions.\nMutation profile analysis: We constructed a binary mutation matrix (presence/absence) for all protein-\ncoding genes across all patients. Principal component analysis (PCA) was performed separately for\nmelanoma and NSCLC. To quantify cohort-specific clustering, we calculated silhouette scores for the top\n10 principal components using cohort labels as the clustering variable. Silhouette scores range from -1 to\n1, where values near 0 indicate no meaningful clustering structure, and negative values suggest that cohort\nassignment does not reflect natural groupings in the data, both supporting the absence of cohort-driven\nclustering.\nPerformance evaluation across cohorts: To assess whether batch effects impacted model performance,\nwe conducted within-cohort validation by calculating AUROC for each score separately for each cohort. We\nthen performed 5-fold cross-validation incorporating cohort as a categorical covariate in a logistic regression\nframework.\nFig. R1. Batch effect evaluation across ICI cohorts.\n(A) Pairwise Mann-Whitney U test comparing TMB distributions across cohorts within the same cancer type and\nresponse category for melanoma (left) and NSCLC (right). Heatmaps show FDR-corrected p-values. (new Fig.\nS7a)\n(B) Distribution of top 2 principal components of binary mutation profiles for NSCLC (top) and melanoma (bottom).\nEach point represents a patient, colored by cohort. (new Fig. S7b)\n(C) Cohort-specific AUROC performance comparing TMB, CSiN, ioTNL, NP-LandscapeSum, NP-LandscapeCCF\nand NP-LandscapeClone within individual cohorts. (new Fig. S11b)\n(D) AUROC performance using 5-fold cross-validation repeated 100 times with cohort incorporated as a covariate\nfor melanoma (left) and NSCLC (right). (new Fig. S11c)\n2. The training set that was created was created in mind of finding differences in T cells’ ability to\ndifferentiate peptides with only 1 amino acid difference. How many unique T cells were included per\npMHC triplet? Are these T cells public clones or private clones of patients? My concern about this\nmethod is whether this can be generalized or is it only accounting for specific T cells that are private\na few people?\nThis is an important point regarding model generalizability. To predict similarity in terms of T-cell\nrecognition between wild-type and mutant peptides, we identified cross-reactive peptide pairs from\npMHC-TCR binding datasets. In our study, cross-reactive peptides are defined as peptides presented\nby the same MHC allele and recognized by the same TCR CDR3 sequence. This definition captures\npeptides that are truly cross-reactive and solely due to the change in amino acid composition. Under\nthis strict criterion, we curated a dataset of 1,153 MHC-I and 261 MHC-II cross-reactive peptide pairs.\nAlthough our model does not explicitly incorporate CDR3 sequences, analyzing CDR3 coverage\nprovides important insights into generalizability. We calculated the number of distinct CDR3\nsequences that bind to each cross-reactive peptide pair and found that approximately 50% of\npairs are recognized by multiple CDR3 sequences (Fig. R2A). Since our model focuses on single\namino acid changes, we further examined CDR3 support for peptide pairs sharing the same\namino acid substitution. We found that more than 50% of these substitutions are supported by\nmultiple independent CDR3 sequences (Fig. R2B). This multi-CDR3 support suggests that cross-\nreactivity is driven by generalizable peptide features rather than specific CDR3 sequences. Our model\ncaptures these intrinsic peptide properties by learning from diverse TCR-peptide binding patterns.\nFurthermore, the model performed well at predicting T-cell reactivity in the NCI dataset which is not\nconstrained to the subset of pMHCs or CDR3s on which NeoPrecis was trained and fine-tuned. In\nour response to the next comment, we further address the generalizability across unobserved HLA\nalleles.\nWe acknowledge that due to the vast diversity of the available TCR repertoire (2.5 × 107 - 3.8 × 108),\nour dataset captures only a limited fraction of possible CDR3 sequences. We have added this\nanalysis to the Results section and note that expanding the dataset with additional pMHC-TCR\nbinding data would further strengthen the model's generalizability.\nManuscript updates:\nUpdated Figures: Fig. R2 has been added as Fig. S3a\nUpdated Results: Line 190–196\nTo assess model generalizability, we evaluated TCR diversity in the training data and performance on\nunobserved MHC alleles. For TCR diversity, we analyzed CDR3 coverage in the cross-reactive peptide\ntriplet dataset and found that approximately 50% of cross-reactive peptide pairs are recognized by multiple\ndistinct CDR3 sequences, with more than 50% of amino acid substitutions supported by multiple\nindependent CDR3s (Fig. S3a). This multi-CDR3 support indicates that the training data captures\ngeneralizable cross-reactivity patterns rather than isolated clone-specific interactions.\nFig. R2. Cumulative distribution of the number of distinct CDR3 sequences supporting (A) each cross-reactive\npeptide pair and (B) each amino acid substitution in the triplet dataset. (new Fig. S3a)\n3. Training data size for HLA class II is small. How many alleles were you able to train for? How did you\naccount for overfit? Do you see any skew towards specific alleles for either HLA class I or II?\nWe appreciate this important question about model generalizability across HLA alleles. Due to the\nstrict criteria for the triplet dataset (derived from TCR-pMHC binding data), allele coverage is indeed\nlimited, with only 15 MHC-I and 7 MHC-II alleles represented. However, the fine-tuning dataset\n(immunogenicity data from CEDAR) provides broader coverage with 88 MHC-I and 23 MHC-II alleles.\nTo address the challenge of limited allele diversity, we represent MHC alleles using pMHC binding\nmotifs rather than allele identity, which enables the model to capture a wider range of MHC alleles\nand generalize to unobserved alleles by capturing shared peptide-binding preferences.\nWe evaluated model generalizability across MHC alleles in two complementary ways: performance\non high-frequency versus low-frequency alleles, and performance on unobserved alleles. Both the\ntriplet and CEDAR datasets exhibit allele frequency skew reflecting population-level biases (Fig.\nR3A-B). For instance, HLA-A*02:01, which is prevalent in European/Caucasian populations, is\nsimilarly dominant in our datasets. To assess whether this skew leads to performance bias, we\ncompared model performance on the top 5 most frequent alleles (high-frequency group) versus all\nother alleles (low-frequency group) in the CEDAR test set. This analysis was restricted to MHC-I due\nto sample size limitations for MHC-II. NeoPrecis-Immuno, PRIME, and ICERFIRE outperformed\nNetMHCpan and showed consistent performance across both groups, whereas DeepNeo exhibited\nsubstantially worse performance in the low-frequency allele group (Fig. R3C). These results\ndemonstrate that our motif-based allele representation effectively mitigates potential overfitting to\nover-represented alleles.\nTo directly evaluate generalizability to unobserved alleles, we stratified mutations in the NCI dataset\nby the number of unobserved alleles—those absent from training data. The CEDAR test set contained\nonly 2 unobserved alleles each for MHC-I and MHC-II (Fig. R4A), necessitating the use of NCI for\nthis evaluation. In NCI, 28% of mutations had unobserved MHC-I alleles, while nearly all mutations\nhad at least one unobserved MHC-II allele (Fig. R4B). Although performance decreased in the\nunobserved allele group, NeoPrecis-Immuno outperformed other predictors for MHC-I and performed\ncomparably for MHC-II (Fig. R4C). We note that alleles unobserved in our training data may have\nbeen observed in other predictors' training sets, potentially contributing to their performance in these\ncases. This result shows that our model generalizes effectively across unobserved MHC alleles.\nManuscript updates:\nUpdated Figures: Fig. R4 has been added as Fig. S3b-d\nUpdated Results: Line 196–206\nFor MHC allele generalization, we stratified mutations in the NCI dataset by the number of unobserved\nalleles—those absent from training data. The CEDAR test set contained only 2 unobserved alleles each for\nMHC-I and MHC-II (Fig. S3b), necessitating the use of NCI for this evaluation. In NCI, 28% of mutations\nhad unobserved MHC-I alleles, while nearly all mutations had at least one unobserved MHC-II allele (Fig.\nS3c). Although performance decreased in the unobserved allele group, NeoPrecis-Immuno outperformed\nother predictors for MHC-I and performed comparably for MHC-II (Fig. S3d). We note that alleles\nunobserved in our training data may have been observed in other predictors’ training sets, potentially\ncontributing to their performance in these cases. These results demonstrate that NeoPrecis-Immuno\ngeneralizes effectively across unobserved MHC alleles and that the training data represents diverse TCR\nrecognition patterns, though limited training diversity remains a constraint.\nFig. R3. Performance of high-frequency versus low-frequency alleles.\n(A) Proportion distribution of the top 10 frequent alleles in the triplet dataset.\n(B) Proportion distribution of the top 10 frequent alleles in the CEDAR training set.\n(C) Performance (AUROC and AUPRC) of immunogenicity predictors stratified by allele frequency. High-frequency\nalleles are defined as the top 5 alleles; low-frequency alleles comprise all others. Performance metrics were\ncalculated using bootstrapping (1,000 iterations).\nFig. R4. Performance evaluation of unobserved alleles. (new Fig. S3b-d)\n(A) Venn diagram showing allele overlap between the cross-reactive triplet dataset, CEDAR training set, and\nCEDAR testing set for MHC-I (top) and MHC-II (bottom).\n(B) Cumulative distribution of the number of unobserved alleles per mutation in the NCI dataset for MHC-I (left) and\nMHC-II (right). Unobserved alleles are those absent from the NeoPrecis-Immuno training data.\n(C) Performance (AUROC) of immunogenicity predictors stratified by unobserved allele groups for MHC-I (left) and\nMHC-II (right). Mutations were grouped by the number of unobserved alleles (0 vs. > 0 for MHC-I; ≤ 5 vs. > 5\nfor MHC-II). AUROC distributions were calculated using bootstrapping (1,000 iterations). In the boxplots, the\ncenter line represents the median, boxes indicate the interquartile range (IQR; 25th–75th percentile), and\nwhiskers extend to the most extreme data points within 1.5× the IQR. Outliers are depicted as individual circles.\nMinor comments:\n1. Please rearrange locations of figures 3B, 3C, 3D to make them read from left to right and top to\nbottom for easier readability\nWe thank the reviewer for this suggestion. We have rearranged the layout of Fig. 3b-d to follow a\nleft-to-right, top-to-bottom reading order, which improves the clarity and flow of the figure.\nReviewer #2 (Remarks to the Author):\nThe authors describe a computational framework that predicts the immunogenicity of neoantigens, with\nthe option to include tumor clonality as an additional feature to enhance prediction of response to ICI\ntherapy. The authors show on-par performance of the immunogenicity model when compared with\npreviously published studies using similar approaches. However, the improvement in performance is not\nstrikingly different and could be highly dependent on the test dataset. I believe the novelty of the tool\ncould be more focused on the incorporation of clonality information, and how this additional information\ncan be leveraged to improve predictions across more tumor cohorts. Please see detailed comments\nbelow:\nWe thank this reviewer for the helpful suggestions which have improved the clarity of our work. We have\naddressed all concerns about model generalizability, including MHC allele and peptide length bias,\nthrough comprehensive analyses using unobserved allele and length stratification. Additionally, we have\nemphasized the novel contribution of clonality integration by demonstrating the coordination of MHC-I\nand MHC-II pathways in tumor subclonal immunity, with significant performance improvement when using\nmultiplicative integration of MHC-I and MHC-II predictions at the subclonal level.\nMajor comments:\n1. The authors highlight the use of allele-specific MHC binding motifs as well as position-weighted\npooling as novel features for their immunogenicity prediction model. How were the MHC binding\nmotifs and positional weights derived or calculated? Given that motifs are generally derived from\npeptides that bind or are presented by the respective allele, how does the model handle alleles with\nno training data and thus no clear binding motif?\nWe thank the reviewer for this important question regarding MHC binding motifs and position weights.\nWe acknowledge that our original description of NeoPrecis-Immuno lacked sufficient detail, and we\nhave substantially revised the Methods section to provide a comprehensive explanation of how these\ncomponents are derived and implemented.\nCalculation of MHC binding motifs and positional weights:\nInput preparation: pMHC complexes are first predicted using an MHC-binding predictor, and the 9-\nmer core binding region is identified as the input sequence. For each neoantigen, we obtain paired\nwild-type (WT) and mutant (MT) peptide sequences represented as one-hot encoded matrices:\n𝑃𝑊𝑇,𝑃𝑀𝑇 ∈ℝ𝐿×20, where L=9 is the core binding region length and each row is a one-hot vector over\n20 amino acids.\nMHC binding motif representation: To represent MHC alleles in a generalizable manner, we\nconstructed binding motifs from NetMHCpan predictions. We randomly sampled 100,000 peptides\nfrom the human peptidome (9-mers for MHC-I, 15-mers for MHC-II) and predicted binding affinity for\neach available allele. Peptides passing the binding rank threshold (2 for MHC-I, 10 for MHC-II) were\nretained, and their 9-mer binding cores were used to compute position-specific amino acid\nfrequencies. The resulting binding motif 𝑀 ∈ ℝ𝐿×20 serves as the MHC allele representation, where\neach position contains the frequency distribution over 20 amino acids.\nAmino acid embedding: One-hot encoded amino acids are first transformed using the BLOSUM62\nsubstitution matrix, then projected into a latent space via a single-layer neural network:\n𝐸𝑊𝑇 = 𝑓 (𝑃𝑊𝑇)∈ ℝ𝐿×𝑑𝑎\n𝑒𝑚𝑏𝑒𝑑\n𝐸𝑀𝑇 = 𝑓 (𝑃𝑀𝑇)∈ ℝ𝐿×𝑑𝑎\n𝑒𝑚𝑏𝑒𝑑\n𝐸𝑚𝑜𝑡𝑖𝑓 = 𝑓 (𝑀) ∈ℝ𝐿×𝑑𝑎\n𝑒𝑚𝑏𝑒𝑑\nwhere d =2 is the amino acid embedding dimension.\na\nMotif enrichment: To capture the complex interactions between peptide sequences and MHC alleles,\nwe introduce a motif-dependent convolutional layer. The model learns K=8 convolutional kernels\n𝑊𝑐𝑜𝑛𝑣 ∈ ℝ𝑑𝑎×𝑑𝑎, which are modulated by the MHC motif embedding:\n𝑘\n𝐶 = 𝐸𝑚𝑜𝑡𝑖𝑓𝑊𝑐𝑜𝑛𝑣 ∈ ℝ𝐿×𝑑𝑎,𝑘 = 1,…,𝐾\n𝑘 𝑘\nThese motif-enriched kernels are then applied position-wise to peptide embeddings:\n𝐻𝑊𝑇 = 𝐶 ⨀𝐸𝑊𝑇 ∈ℝ𝐿\n𝑘 𝑘\n𝐻𝑀𝑇 = 𝐶 ⨀𝐸𝑀𝑇 ∈ ℝ𝐿\n𝑘 𝑘\nwhere ⨀ denotes element-wise multiplication followed by summation across the embedding\ndimension. The K kernel responses are concatenated and projected to the peptide embedding\ndimension:\n𝐸𝑊𝑇 = 𝑊𝑝𝑟𝑜𝑗[𝐻𝑊𝑇,…,𝐻𝑊𝑇] ∈ℝ𝐿×𝑑𝑝\n𝑒𝑛𝑟𝑖𝑐ℎ𝑒𝑑 1 𝐾\n𝐸𝑀𝑇 = 𝑊𝑝𝑟𝑜𝑗[𝐻𝑀𝑇,…,𝐻𝑀𝑇]∈ ℝ𝐿×𝑑𝑝\n𝑒𝑛𝑟𝑖𝑐ℎ𝑒𝑑 1 𝐾\nwhere 𝑊𝑝𝑟𝑜𝑗 ∈ℝ𝑑𝑝×𝐾 is a learnable projection matrix and d =2 is the peptide embedding dimension.\np\nPosition-weighted peptide pooling: To aggregate position-enriched embeddings into fixed-length\npeptide representations, we employ position-weighted pooling. Each of the 9 core binding positions\nhas a learnable weight w (separate parameters for MHC-I and MHC-II):\ni\n𝐿\n𝑣𝑊𝑇 = ∑𝑤 𝐸𝑊𝑇 ∈ ℝ𝑑𝑝\n𝑖 𝑒𝑛𝑟𝑖𝑐ℎ𝑒𝑑,𝑖\n𝑖=1\n𝐿\n𝑣𝑀𝑇 = ∑𝑤 𝐸𝑀𝑇 ∈ ℝ𝑑𝑝\n𝑖 𝑒𝑛𝑟𝑖𝑐ℎ𝑒𝑑,𝑖\n𝑖=1\nwhere 𝐸 denotes the embedding at position i. This design allows the model to emphasize\n𝑒𝑛𝑟𝑖𝑐ℎ𝑒𝑑,𝑖\nstructurally important positions in the MHC binding groove.\nGeometric representation: To capture the distinction between wild-type and mutant peptides in a\nway that reflects TCR cross-reactivity, we construct a geometric representation in the peptide\nembedding space. The geometric transformation consists of three components:\n𝑣𝑜𝑟𝑖𝑔𝑖𝑛 = 𝑣𝑊𝑇 ∈ ℝ𝑑𝑝\n𝑣𝑀𝑇−𝑣𝑊𝑇\n𝑣𝑑𝑖𝑟𝑒𝑐𝑡𝑖𝑜𝑛 = ∈ ℝ𝑑𝑝\n‖𝑣𝑀𝑇−𝑣𝑊𝑇‖ +𝜀\n2\nd = ‖𝑣𝑀𝑇 −𝑣𝑊𝑇‖ ∈ ℝ\n2\nwhere vorigin represents the wild-type peptide position, vdirection is the normalized direction of mutation\nin the embedding space, d is the Euclidean distance between wild-type and mutant embeddings,\nand ε is a small constant to prevent division by zero.\nCross-reactivity distance (CRD): To account for sequence context, we compute a scaling factor\nusing a two-layer neural network:\n𝑠 = 𝑓 ([𝑣𝑜𝑟𝑖𝑔𝑖𝑛;𝑣𝑑𝑖𝑟𝑒𝑐𝑡𝑖𝑜𝑛])\n𝑠𝑐𝑎𝑙𝑒\nThe CRD is then computed as:\n𝐶𝑅𝐷 = 𝑑∙𝜎(𝑠)\nwhere σ is the sigmoid function. Higher CRD values indicate lower predicted TCR cross-reactivity\nbetween wild-type and mutant peptides.\nImmunogenicity prediction: The final immunogenicity prediction integrates CRD with MHC binding\nprediction:\n𝑦 = 𝜎(𝑊 [𝐶𝑅𝐷;𝑀𝐻𝐶 ]+𝑏 )\n𝑖𝑚𝑚𝑢𝑛𝑜 𝑜𝑢𝑡 𝑏𝑖𝑛𝑑𝑖𝑛𝑔 𝑜𝑢𝑡\nThis produces an immunogenicity score between 0 and 1.\nModel training: The NeoPrecis-Immuno model was trained in two stages. In the first stage, it was\ntrained on the cross-reactive peptide triplet dataset using contrastive learning with triplet loss to\noptimize CRD:\nℒ = 𝑚𝑎𝑥(0,𝐶𝑅𝐷(𝑎,𝑝)−𝐶𝑅𝐷(𝑎,𝑛)+𝑚𝑎𝑟𝑔𝑖𝑛)\n𝑡𝑟𝑖𝑝𝑙𝑒𝑡\nwhere a represents an anchor peptide, p is a cross-reactive peptide, n is a non-cross-reactive\npeptide, and the margin was set to 1. In the second stage, the model was fine-tuned on the CEDAR\nimmunogenicity dataset using binary cross-entropy loss:\n1\nℒ = − [𝑦 log(𝑦̂ )+(1−𝑦 )log(1−𝑦̂ )]\n𝐵𝐶𝐸 𝑁 𝑖 𝑖 𝑖 𝑖\nwhere 𝑦 is the true immunogenicity label and 𝑦̂ is the predicted score. During fine-tuning, all layers\n𝑖 𝑖\nwere frozen, except for the scaling factor and classifier, allowing the model to focus on adjusting\nimmunogenicity predictions. Training was conducted for 100 epochs with a learning rate of 0.005,\nusing the AdamW optimizer.\nUnobserved alleles:\nThe reviewer also raises a critical consideration regarding alleles absent from the training data. While\nour motif-based representation is designed to generalize across alleles, evaluating performance on\nunobserved alleles is essential to validate this approach.\nAtkins et al. demonstrated that NetMHCpan maintains strong performance on alleles not included in\nits training set, which supports the reliability of binding motifs derived from NetMHCpan predictions.\nBuilding on this foundation, we directly evaluated NeoPrecis-Immuno performance on unobserved\nalleles. The CEDAR test set contained only 2 unobserved alleles each for MHC-I and MHC-II (Fig.\nR4A), necessitating the use of the NCI dataset for this evaluation.\nSince the NCI dataset is a mutation-level dataset where each mutation can potentially bind to any of\nan individual's MHC alleles, we calculated the proportion of unobserved alleles per mutation. We\nfound that most MHC-I alleles were represented in the training data, with only 28% of mutations\ninvolving unobserved MHC-I alleles, whereas nearly all mutations had at least one unobserved MHC-\nII allele (Fig. R4B). We then stratified mutations by the number of unobserved alleles and compared\nmodel performance. Although performance decreased in the unobserved allele group, NeoPrecis-\nImmuno outperformed other predictors for MHC-I and performed comparably for MHC-II (Fig. R4C).\nWe note that alleles unobserved in our training data may have been observed in other predictors'\ntraining sets, potentially contributing to their performance in these cases. These results demonstrate\nthat our motif-based representation successfully generalizes to unobserved alleles despite limited\ntraining data for these alleles.\nManuscript updates:\nUpdated Figures: Fig. R4 has been added as Fig. S3b-d\nUpdated Results: Line 196–206\nFor MHC allele generalization, we stratified mutations in the NCI dataset by the number of unobserved\nalleles—those absent from training data. The CEDAR test set contained only 2 unobserved alleles each for\nMHC-I and MHC-II (Fig. S3b), necessitating the use of NCI for this evaluation. In NCI, 28% of mutations\nhad unobserved MHC-I alleles, while nearly all mutations had at least one unobserved MHC-II allele (Fig.\nS3c). Although performance decreased in the unobserved allele group, NeoPrecis-Immuno outperformed\nother predictors for MHC-I and performed comparably for MHC-II (Fig. S3d). We note that alleles\nunobserved in our training data may have been observed in other predictors’ training sets, potentially\ncontributing to their performance in these cases. These results demonstrate that NeoPrecis-Immuno\ngeneralizes effectively across unobserved MHC alleles and that the training data represents diverse TCR\nrecognition patterns, though limited training diversity remains a constraint.\nUpdated Methods: The details of NeoPrecis model architecture have been add to Line 570–646\nReference:\nAtkins TK, Solanki A, Vasmatzis G, Cornette J, Riedel M. Evaluating NetMHCpan performance on non-\nEuropean HLA alleles not present in training data. Front Immunol. 2024 Jan 16;14:1288105. doi:\n10.3389/fimmu.2023.1288105. PMID: 38292493; PMCID: PMC10825027.\nFig. R4. Performance evaluation of unobserved alleles. (new Fig. S3b-d)\n(A) Venn diagram showing allele overlap between the cross-reactive triplet dataset, CEDAR training set, and\nCEDAR testing set for MHC-I (top) and MHC-II (bottom).\n(B) Cumulative distribution of the number of unobserved alleles per mutation in the NCI dataset for MHC-I (left) and\nMHC-II (right). Unobserved alleles are those absent from the NeoPrecis-Immuno training data.\n(C) Performance (AUROC) of immunogenicity predictors stratified by unobserved allele groups for MHC-I (left) and\nMHC-II (right). Mutations were grouped by the number of unobserved alleles (0 vs. > 0 for MHC-I; ≤ 5 vs. > 5\nfor MHC-II). AUROC distributions were calculated using bootstrapping (1,000 iterations). In the boxplots, the\ncenter line represents the median, boxes indicate the interquartile range (IQR; 25th–75th percentile), and\nwhiskers extend to the most extreme data points within 1.5× the IQR. Outliers are depicted as individual circles.\n2. Distance based on BLOSUM62 was used throughout the manuscript for multiple comparisons (e.g.,\ncomparing to the geometric distance calculated by NeoPrecis-Immuno). Could the authors show\nsimilar comparisons using PMBEC, or further explain why PMBEC was not considered a more\nbiologically relevant distance matrix?\nWe thank the reviewer for this insightful suggestion to compare our model with PMBEC, an amino\nacid substitution matrix specifically derived from peptide-MHC binding data. We have now\nincorporated PMBEC comparisons into our analysis.\nComparison of amino acid representations:\nWe obtained the PMBEC encoding from Kim et al. and compared three amino acid representation\nmethods: BLOSUM62, PMBEC, and NeoPrecis-Immuno embeddings. For BLOSUM62 and PMBEC,\nwe applied PCA to extract the first two principal components as two-dimensional embeddings for\nvisualization (Fig. R5A). The three representation methods yield distinctly different embedding\nstructures, reflecting their different derivation approaches and biological contexts. For example, in\nNeoPrecis-Immuno, histidine (H) was positioned closer to tyrosine (Y) and phenylalanine (F), likely\ndue to their shared aromatic rings, while tryptophan (W) was positioned further away, possibly\nreflecting its larger size. These differences between NeoPrecis and PMBEC suggests distinct amino\nacid interaction patterns in TCR recognition versus MHC presentation.\nPerformance comparison on immunogenicity prediction:\nWe evaluated the predictive power of amino acid distances derived from each method for\nimmunogenicity prediction, following the approach in Fig. 3a (Fig. R5B). We tested two distance\ncalculation strategies: (1) Euclidean distance in PCA-reduced space (top 2 PCs) to obtain continuous\ndistance distributions (BLOSUMDist, PMBECDist), and (2) direct substitution scores from the\nmatrices (BLOSUMDist(sub), PMBECDist(sub)). The substitution distance derived from NeoPrecis-\nImmuno (SubDist) outperformed all baseline metrics regardless of distance calculation method.\nManuscript updates:\nIn the main manuscript, we report only the PC-derived distance metrics (BLOSUMDist, PMBECDist;\nFig. 3a) to maintain continuous distance distributions, which is particularly important for Fig. S2a. We\nhave added PMBEC as a baseline comparison and updated the Methods section to clarify the amino\nacid distance calculation approach.\nUpdated Figures: Fig. 3a (component analysis), Fig. S2a (distribution of PMBEC distance), Fig. S4a\n(Correlation between components)\nUpdated Results: Line 211–212\nBLOSUM62 (BLOSUMDist) and PMBEC (PMBECDist) distances were included as baselines for\ncomparison.\nUpdated Methods: Line 680–686\nFor baseline comparison, we used BLOSUM62 and PMBEC substitution matrices. BLOSUM62 is an amino\nacid substitution matrix derived from aligned protein sequences, while PMBEC is specifically derived from\npeptide-MHC binding data. Since substitution matrices provide discrete scores for amino acid pairs, we\napplied PCA to each matrix to extract the first two principal components as two-dimensional amino acid\nembeddings, enabling continuous distance calculations. We then computed distances as the Euclidean\ndistance between wild-type and mutated amino acids in this reduced space (BLOSUM62Dist, PMBECDist).\nReference:\nKim Y, Sidney J, Pinilla C, Sette A, Peters B. Derivation of an amino acid similarity matrix for peptide: MHC\nbinding and its application as a Bayesian prior. BMC Bioinformatics. 2009 Nov 30;10:394. doi: 10.1186/1471-\n2105-10-394. PMID: 19948066; PMCID: PMC2790471.\nFig. R5. Comparison of amino acid embeddings and distance metrics.\n(A) Two-dimensional amino acid embeddings derived from BLOSUM62 (left), PMBEC (center), and NeoPrecis\n(right). For BLOSUM62 and PMBEC, the top 2 principal components from PCA are shown. NeoPrecis\nembeddings are learned during model training.\n(B) AUROC comparison of amino acid distance metrics and incrementally constructed NeoPrecis-Immuno\ncomponents for MHC-I (left) and MHC-II (right) predictions on the NCI dataset. Distance metrics include:\nBLOSUM_Dist (Euclidean distance in top 2 PC space), BLOSUM_Dist(sub) (negative substitution score),\nPMBEC_Dist (Euclidean distance in top 2 PC space), PMBEC_Dist(sub) (negative substitution score),\nSubDist (NeoPrecis substitution distance), SubPosDist (SubDist with position weighting), GeoDist (SubPosDist with\nmotif enrichment), and CRD (GeoDist with sigmoid scaling). (new Fig. 3a)\n3. For NeoPrecis-Immuno’s evaluation on the CEDAR dataset (internal test set), the positive-to-negative\nratio is roughly 1:1 (line 165). This ratio seems quite different from other benchmark datasets (e.g.,\ntest sets used in ICERFIRE and PRIME). Does this ratio reflect the true proportion of peptides that\ncan be recognized by T cells? Additionally, for the external test set, the authors mention aggregation\nmethods across multiple peptides for determining mutation site immunogenicity. Have the authors\ntried and compared multiple aggregation strategies, and are performance results dependent on the\naggregation method? Do the authors mask positive peptides found in training prior to aggregation\nand evaluation? Could the authors provide and comment on model performance comparisons at an\nallele+length-specific level?\nWe appreciate these important concerns about real-world applicability of neoantigen prediction\nmodels. We have conducted additional analyses and revised the manuscript to address the positive-\nto-negative ratio, aggregation strategies, and allele- and length-specific performance.\nPositive-to-negative ratio:\nIn this study, we aimed to predict T-cell recognition of MHC-presented peptides. However, the true\npositive rate is difficult to estimate and varies considerably across studies due to differences in\nexperimental methods and selection criteria.\nThe experimental workflow significantly impacts observed positive rates. Studies following a \"top-\ndown\" approach—first selecting MHC-binding peptides, then testing for T-cell responses—typically\nobserve positive rates of 20-40%. In the CEDAR dataset, the overall positive rate using\nmultimer/tetramer assays is 22.2% for MHC-I (493/2,223), while Alban et al. reported 39.0% (196/502)\nin non-small cell lung cancer using a similar top-down strategy. In contrast, \"bottom-up\" approaches\nthat perform T-cell assays without MHC binding pre-selection yield much lower rates. The NCI dataset\nexemplifies this: ELISPOT screening without MHC filtering yielded only 0.6% positive (41/7,384),\nincreasing to 3.3% (36/1,088) when restricted to predicted MHC binders. Our model follows the\nbiologically sound top-down principle by integrating MHC binding with T-cell cross-reactivity metrics.\nWe acknowledge that the 50% positive rate in our CEDAR test set does not reflect realistic clinical\nscenarios, resulting from filtering for predicted MHC binders and selecting neoantigens absent from\nother predictors' training sets. For fair comparison with existing immunogenicity predictors, we\nexamined the test set compositions used in their publications: PRIME used 25.6% positive\n(1,282/4,958) and ICERFIRE used 20.8% positive (631/3,033). To match these benchmarks, we\nrebalanced samples through bootstrapping to adjust our positive rate to 25%. The resulting AUROC\nand AUPRC rankings remained consistent across predictors (Fig. R6A), demonstrating that our\nmodel's performance is robust to class imbalance.\nAggregation strategies:\nTo predict mutation-level immunogenicity from allele-specific scores, we evaluated multiple\naggregation strategies that account for individual MHC genotypes. Each strategy involves two key\ndecisions: whether to apply allele masking and which aggregation method to use.\nAllele masking: Masking excludes alleles predicted not to bind the target peptide, based on\nNetMHCpan %rank thresholds (≤ 2 for MHC-I, ≤ 10 for MHC-II). Let A denote all MHC alleles for\nan individual and 𝐴 ⊆ 𝐴 denote binding alleles. Without masking, all alleles contribute to\naggregation; with masking, only A is considered.\nAggregation methods: For allele-specific scores s , we evaluated five aggregation methods.\n- Mean: 𝑠 = ∑ 𝑠\n𝑚𝑢𝑡 𝑎∈𝐴 𝑎\n|𝐴|\n- Harmonic mean: 𝑠 = |𝐴|\n𝑚𝑢𝑡 ∑ 𝑎∈𝐴(𝑠𝑎+𝜀)−1\n- Weighted average: 𝑠 = ∑ 𝑎∈𝐴𝑤𝑎𝑠𝑎,\nwhere w is binding affinity from NetMHCpan\n𝑚𝑢𝑡 ∑ 𝑎∈𝐴𝑤𝑎 a\n- Best binding: 𝑠 = 𝑠 , where 𝑎∗ = 𝑎𝑟𝑔𝑚𝑖𝑛 𝑟𝑎𝑛𝑘(𝑎)\n𝑚𝑢𝑡 𝑎∗ 𝑎∈𝐴\n- Maximum: 𝑠 = 𝑚𝑎𝑥 𝑠\n𝑚𝑢𝑡 𝑎∈𝐴 𝑎\nWe compared all ten combinations (with/without masking × five aggregation methods) on the NCI\ndataset and found that masked maximum approach performed most consistently across both MHC\nclasses and different immunogenicity metrics (Fig. R6B; Fig. S5b).\nAllele- and length-specific performance\nOur model generalizes across alleles and peptide lengths by leveraging NetMHCpan predictions. We\nrepresent alleles using binding motifs and standardize inputs to 9-mer core binding regions regardless\nof original peptide length. However, to assess potential biases across different alleles and lengths,\nwe conducted peptide-level performance comparisons on the CEDAR test set.\nDue to limited sample sizes, we analyzed alleles and lengths separately. We identified alleles and\nlengths with at least 40 peptides, then used bootstrapping to compute performance metrics across\nmodels. DeepNeo was excluded from this analysis as it only accepts 9-mer inputs.\nCompared to the NetMHCpan baseline, NeoPrecis-Immuno demonstrated the most consistent\nperformance, outperforming the baseline in 4 of 5 alleles and both lengths (Fig. R6C). In contrast,\nPRIME and ICERFIRE showed substantial variation across alleles, with PRIME outperforming the\nbaseline in only 2 of 5 alleles and ICERFIRE in only 1 allele. These results indicate that NeoPrecis-\nImmuno exhibits less allele-specific bias and maintains robust performance across diverse MHC\ncontexts.\nManuscript updates:\nUpdated Figures: Fig. S2d (performance of rebalanced testing set), Fig. S5b (aggregation analysis)\nUpdated Results: Line 177–180\nWe also rebalanced samples through bootstrapping to achieve a 1:3 positive-to-negative ratio, matching\nthe test set balance used in PRIME and ICERFIRE benchmarks and better reflecting realistic\nimmunogenicity rates, which yielded similar performance rankings (Fig. S2d).\nUpdated Methods: The above aggregation method has been updated to Line 748–764\nReference:\nAlban TJ, Riaz N, Parthasarathy P, Makarov V, Kendall S, Yoo SK, Shah R, Weinhold N, Srivastava R, Ma X,\nKrishna C, Mok JY, van Esch WJE, Garon E, Akerley W, Creelan B, Aanur N, Chowell D, Geese WJ, Rizvi NA,\nChan TA. Neoantigen immunogenicity landscapes and evolution of tumor ecosystems during immunotherapy\nwith nivolumab. Nat Med. 2024 Nov;30(11):3209-3222. doi: 10.1038/s41591-024-03240-y. Epub 2024 Sep 30.\nPMID: 39349627; PMCID: PMC12066197.\nFig. R6. Evaluation of NeoPrecis-Immuno under class rebalance, aggregation approach, and allele or length\nstratification.\n(A) Performance (AUROC and AUPRC) for each predictor on rebalanced MHC-I samples (n = 438; positive-to-\nnegative ratio adjusted to 1:3 through bootstrapping). Performance distributions were calculated using\nbootstrapping (1,000 iterations). (new Fig. S2d)\n(B) The heatmap displays the negative log p-values from the association tests (Mann–Whitney U test) between\neach approach and T-cell activity across various TCR recognition metrics. Annotations indicate the rank within\neach MHC class and for each metric. (new Fig. S5b)\n(C) AUROC of each predictor stratified by allele (left) or peptide length (right) for MHC-I samples with ≥40 data\npoints per group. Performance was evaluated using bootstrapping (1,000 iterations).\n4. Given the relatively small dataset sizes, could the authors clarify the number of alleles and peptide\nlengths present in the training/testing data? How reliable is the model when generalizing to alleles\nnot seen during training (do such cases exist in the test set)?\nWe thank the reviewer for this important question about dataset composition and model\ngeneralizability to unobserved alleles.\nDataset composition:\nDue to strict criteria for identifying cross-reactive peptide pairs in the triplet training dataset, allele\ncoverage is indeed limited, with only 15 MHC-I and 7 MHC-II alleles represented. However, the fine-\ntuning dataset (immunogenicity data from CEDAR) provides substantially broader coverage with 88\nMHC-I and 23 MHC-II alleles (Fig. R4A in response to Reviewer 2 Comment 1).\nFor the peptide length, NeoPrecis uses the binding core region derived from NetMHCpan as input,\nwhere the binding core length is consistently 9-mer regardless of the original peptide length. To\nassess potential length bias, we examined the distribution of original peptide lengths in our datasets.\nIn the CEDAR training set, 9-mer peptides predominate for MHC-I alleles, while 15-mer peptides are\nmost common for MHC-II alleles (Fig. R7). Since DeepNeo only accepts 9-mer peptides, we focused\non 9-mer performance in the main figures (Fig. 2c). For a more comprehensive length evaluation, we\nanalyzed both 9-mer and 10-mer peptides in the CEDAR test set as described in our response to the\nprevious comment (Fig. R6C in response to Reviewer 2 Comment 3). Comparison with the\nNetMHCpan baseline revealed no obvious performance bias across these peptide lengths.\nEvaluation on unobserved alleles:\nThe CEDAR test set contains only 2 unobserved alleles each for MHC-I and MHC-II (Fig. R4A in\nresponse to Reviewer 2 Comment 1), limiting our ability to assess generalization within this dataset.\nTo more comprehensively evaluate performance on unobserved alleles, we conducted additional\nanalysis using the independent NCI dataset, as described in our response to Comment 1.\nWe stratified mutations by the number of unobserved alleles and compared model performance\nacross groups. We found that most MHC-I alleles were represented in the training data, with only\n28% of mutations involving unobserved MHC-I alleles, whereas nearly all mutations had at least one\nunobserved MHC-II allele (Fig. R4B in response to Reviewer 2 Comment 1). Although performance\ndecreased in the unobserved allele group, NeoPrecis-Immuno outperformed other predictors for\nMHC-I and performed comparably for MHC-II (Fig. R4C in response to Reviewer 2 Comment 1).\nWe note that alleles unobserved in our training data may have been observed in other predictors'\ntraining sets, potentially contributing to their performance in these cases. This result shows that our\nmodel generalizes effectively across unobserved MHC alleles.\nManuscript updates:\nUpdated Figures: Fig. R4 has been added as Fig. S3b-d\nUpdated Results: Line 196–206\nFor MHC allele generalization, we stratified mutations in the NCI dataset by the number of unobserved\nalleles—those absent from training data. The CEDAR test set contained only 2 unobserved alleles each for\nMHC-I and MHC-II (Fig. S3b), necessitating the use of NCI for this evaluation. In NCI, 28% of mutations\nhad unobserved MHC-I alleles, while nearly all mutations had at least one unobserved MHC-II allele (Fig.\nS3c). Although performance decreased in the unobserved allele group, NeoPrecis-Immuno outperformed\nother predictors for MHC-I and performed comparably for MHC-II (Fig. S3d). We note that alleles\nunobserved in our training data may have been observed in other predictors’ training sets, potentially\ncontributing to their performance in these cases. These results demonstrate that NeoPrecis-Immuno\ngeneralizes effectively across unobserved MHC alleles and that the training data represents diverse TCR\nrecognition patterns, though limited training diversity remains a constraint.\nFig. R7. Distribution of original peptide lengths in the CEDAR training set.\n5. When calculating a tumor-centric immunogenicity score using NeoPrecis-Immuno predictions (line\n298), the authors directly multiply MHC-I and MHC-II predictions for the same mutation. This seems\nto prioritize a single mutation that activates both a class I and class II response, rather than two\nseparate mutations where one activates class I and the other activates class II. Can the authors\nelaborate on why they believe this is the appropriate strategy?\nWe thank the reviewer for this thoughtful question about our strategy for integrating MHC-I and MHC-\nII immunogenicity predictions. We evaluate the MHC integration at both the mutation level and the\nsubclone level.\nMutation level:\nOur multiplication approach prioritizes mutations that activate both CD4+ (MHC-II) and CD8+ (MHC-\nI) T-cell responses. Evidence demonstrates that CD8+ T cells require CD4+ T-cell help at many\nstages of differentiation and maintenance through T-T cooperation based on associative recognition\nof antigens. In further support of what our group believes to be a tenet of protective T cell immunity,\na 2024 paper by Espinosa-Carrasco et al. reported that effective anti-tumor immunity requires a\nfunctional triad of CD4+ T cells, CD8+ T cells, and dendritic cells. Critically, this coordinated response\ndepends on the presentation of both CD4 and CD8 neoantigens within the same antigen presenting\ncell. Under this model, a single mutation capable of activating both MHC-I and MHC-II responses is\nmore likely to elicit coordinated immunity across the tumor cell population. The multiplication strategy\nreflects this biological principle by assigning higher scores to mutations with dual immunogenic\npotential.\nTo determine the optimal integration strategy, we compared three aggregation methods:\nmultiplication (I×II), summation (I+II), and maximum (max(I, II)). All methods show comparable\nperformance in AUROC and AUPRC (Fig. R8A-B) at the mutation level.\nSubclone level:\nWe evaluated the same three MHC-I and MHC-II integration strategies at the subclone level. Clonal\nanalysis was performed using PyClone, which grouped mutations into clusters representing tumor\nsubclones. Within each cluster, NeoPrecis-Immuno scores were summed separately to derive\ncluster-specific immunogenicity for MHC-I and MHC-II, which were then integrated using\nmultiplication, summation, or maximum. The MHC-integrated cluster immunogenicity scores were\nthen averaged across clusters, weighted by their prevalence, to calculate the tumor-centric\nimmunogenicity score.\nMultiplication consistently outperformed other methods in both melanoma and NSCLC (Fig. R8C-D).\nThis superior performance at the subclone level provides additional biological support for the\ncoordination between CD4+ T cell, CD8+ T cell, and dendritic cell. Since tumor cells within the same\nsubclone share mutations, dendritic cells in that subclone can present multiple neoantigens on both\nMHC-I and MHC-II, enabling coordinated dual-class T-cell responses against the same tumor\nsubclone.\nConclusion:\nGiven that multiplication outperformed other approaches at subclone levels, and considering the\nbiological principle that coordinated CD4+ and CD8+ responses enhance anti-tumor immunity, we\nselected multiplication as our integration strategy. This approach prioritizes mutations capable of\nactivating both MHC pathways, which may be particularly important in tumors where such coordinated\nresponses are critical for therapeutic efficacy.\nManuscript updates:\nUpdated Figures: Fig. R8A-B has been added as Fig. S8b; Fig. R8C-D has been added as Fig. S10.\nUpdated Results:\nLine 340–358:\nTo conduct a more comprehensive mutation-level evaluation, we calculated a tumor-centric immunogenicity\nscore, NeoPrecis-LandscapeSum, by summing NeoPrecis-Immuno predictions across all mutations\n(Methods). For each mutation, the MHC-dual NeoPrecis-Immuno score was derived by multiplying the\nMHC-I and MHC-II NeoPrecis-Immuno predictions. This multiplication strategy prioritizes mutations\ncapable of activating both CD4+ and CD8+ T-cell responses. Effective anti-tumor immunity requires\ncoordinated CD4+ T-cell, CD8+ T-cell, and dendritic cell responses, which depends on presenting both\nMHC-I and MHC-II neoantigens within the same antigen-presenting cell43. Given intratumoral heterogeneity,\na single mutation activating both pathways is more likely to elicit coordinated immunity than relying on\nseparate mutations that may not co-occur in the same cell. We compared multiplication (I×II) with\nsummation (I+II) and maximum (max(I, II)) aggregation methods. While all methods showed comparable\nperformance at the mutation level (Fig. S8b), multiplication was selected based on its superior performance\nat the subclone level (see below, Fig. S10). This integrated approach demonstrated superior performance\nover TMB and TNB in both melanoma (AUROC = 0.66) and NSCLC (AUROC = 0.75) (Fig. 5d). Given the\nbetter performance of MHC-dual multiplication, we refer to the tumor-centric immunogenicity score as MHC-\ndual multiplication in subsequent analyses. In survival analyses (Methods), patients in the top half of the\ntumor-centric score showed significantly better outcomes than those in the bottom half (Fig. 5e). NeoPrecis-\nLandscapeSum outperformed TMB in melanoma and showed comparable performance in NSCLC (Fig.\nS8c).\nLine 378–385:\nWe evaluated the same three MHC-I and MHC-II integration strategies (multiplication, summation,\nmaximum) at the subclone level, where subclone-level MHC-I and MHC-II scores were calculated\nseparately then integrated (Methods). Multiplication consistently outperformed other methods in both\nmelanoma and NSCLC (Fig. S10). This superior performance at the subclone level provides additional\nbiological support for the coordination between CD4+ T cell, CD8+ T cell, and dendritic cell. Since tumor\ncells within the same subclone share mutations, dendritic cells in that subclone can present multiple\nneoantigens on both MHC-I and MHC-II, enabling coordinated dual-class T-cell responses against the\nsame tumor subclone.\nUpdated Discussion: Line 479-483\nThese findings align with our results, showing that combining MHC-I and MHC-II predictions improves ICI\nresponse prediction at the subclone level (coordinated immune responses against tumor subclones\npresenting multiple neoantigens on both MHC classes; Fig. S10) and shows comparable performance at\nthe mutation level (single mutations activating dual pathways; Fig. 5b, Fig. S8b).\nUpdated Methods: Line 828-834\nThe NeoPrecis-LandscapeClone method integrated tumor clonality to offer a more nuanced perspective.\nUsing PyClone, mutations were grouped into clusters according to their prevalence in the tumor and\nmutation size. Within each cluster, NeoPrecis-Immuno scores were summed separately to derive cluster-\nspecific immunogenicity for MHC-I and MHC-II, which were then integrated using multiplication (summation\nand maximum were also evaluated for comparison). The MHC-integrated, cluster-centric immunogenicity\nscores were then averaged across clusters, weighted by their prevalence, to calculate the tumor-centric\nimmunogenicity score (Fig. 6a).\nReference:\nEspinosa-Carrasco G, Chiu E, Scrivo A, Zumbo P, Dave A, Betel D, Kang SW, Jang HJ, Hellmann MD, Burt\nBM, Lee HS, Schietinger A. Intratumoral immune triads are required for immunotherapy-mediated elimination\nof solid tumors. Cancer Cell. 2024 Jul 8;42(7):1202-1216.e8. doi: 10.1016/j.ccell.2024.05.025. Epub 2024 Jun\n20. PMID: 38906155; PMCID: PMC11413804.\nFig. R8. Performance evaluation of mutation-level and subclone-level combination of MHC-I and MHC-II\nimmunogenicity.\n(A) AUROC comparison of NeoPrecis-LandscapeSum scores (mutation level) using MHC-I-only, MHC-II-only, or\nMHC-dual (combined MHC-I and MHC-II). For MHC-dual, three aggregation methods were compared:\nmultiplication (I×II), summation (I+II), and maximum (max(I,II)). (new Fig. S8b)\n(B) AUPRC comparison of NeoPrecis-LandscapeSum scores (mutation level). (new Fig. S8b)\n(C) AUROC comparison of NeoPrecis-LandscapeClone scores (subclone level). (new Fig. S10)\n(D) AUPRC comparison of NeoPrecis-LandscapeClone scores (subclone level). (new Fig. S10)\n6. The authors primarily focused on demonstrating the incorporation of tumor clonality information for\ngenerating a tumor-level overall score. However, it would be interesting to see whether incorporating\nclonality at the mutation level (without using DNA VAF as a proxy) could help improve neoantigen\nimmunogenicity predictions.\nWe thank the reviewer for this excellent suggestion to explore mutation-level clonality. We obtained\ncancer cell fraction (CCF) data from Muller et al., for the NCI dataset and conducted additional\nanalyses to assess whether CCF improves mutation-level immunogenicity prediction compared to\nDNA VAF.\nWe first examined the correlation between CCF and DNA VAF. The Pearson correlation was only\n0.236 (Fig. R9A), suggesting that DNA VAF may not be an adequate proxy for mutation clonality,\nlikely due to confounding factors such as copy number alterations and tumor purity.\nWe evaluated CCF in two contexts. First, we replaced DNA VAF with CCF in the NeoPrecis-\nIntegrated model—a logistic regression model incorporating clonality, RNA VAF, RNA expression,\nPHBR-II, and NeoPrecis-Immuno scores as features. Despite the low correlation between DNA VAF\nand CCF, cross-validation performance was comparable between the two models (Fig. R9B),\nsuggesting that clonality information contributes similarly regardless of whether it is measured by\nCCF or DNA VAF.\nSecond, we tested multiplicative integration by combining NeoPrecis-Immuno scores with either CCF\nor DNA VAF, analogous to our tumor-level approach. Neither CCF nor DNA VAF integration improved\nupon NeoPrecis-Immuno alone (Fig. R9C). These findings indicate that mutation-level clonality\nprovides limited additional benefit for immunogenicity prediction when high-quality immunogenicity\nscores are already available, in contrast to tumor-level analyses where clonality effectively prioritizes\nclonal neoantigens.\nFig. R9. Incorporating clonality at the mutation level.\n(A) Correlation between DNA AF and CCF in the NCI dataset. r is Pearson correlation coefficient.\n(B) Comparison of integrated models (A+P+R) using either DNA VAF or CCF as the DNA abundance metric.\nPerformance was evaluated using 4-fold cross-validation repeated 100 times.\n(C) Performance (AUROC and AUPRC) of multiplicative integration combining NeoPrecis-Immuno scores with\neither CCF or DNA AF.\nMinor comments:\n1. When describing the construction of the triplet dataset, it was a bit confusing that the authors describe\nthe positives as non-immunogenic mutated peptides and the negatives as immunogenic (lines 145–\n147), while in Figure S1A, the seed and positive peptides are described as having corresponding\nTCRs. Given that negative peptides were randomly selected, it is unclear why they are described as\nimmunogenic mutated peptides.\nWe appreciate the opportunity to clarify this source of confusion. The underlying principle of our pre-\ntraining approach is to model the distance between wild-type and mutant peptides in terms of T-cell\nrecognition. We leverage cross-reactivity in peptide-TCR interactions to represent peptide similarity:\ncross-reactive peptides (those recognized by the same TCR) are considered similar, while peptides\nnot sharing TCR recognition are considered dissimilar. This learned distance metric is then applied\nto immunogenicity prediction, where mutant peptides similar to their wild-type counterparts (high\ncross-reactivity) are predicted to be less immunogenic, while dissimilar mutant peptides (low cross-\nreactivity) are predicted to be more immunogenic.\nIn the original manuscript, we used \"positive\" to refer to peptides that bind the same TCR as the seed\npeptide, which indeed represents close distance and lower predicted immunogenicity. We recognize\nthis terminology was confusing, as \"positive\" could be misinterpreted as \"immunogenic.\"\nWe have revised the terminology throughout the manuscript for clarity:\n\"TCR-binding triplet\" → \"Cross-reactive peptide triplet\"\n\"Positive peptide\" → \"Cross-reactive peptide\" (shares TCR binding with seed)\n\"Negative peptide\" → \"Non-cross-reactive peptide\" (does not share TCR binding)\nUnder this revised terminology, Figure S1a now clearly shows that the seed and cross-reactive\npeptides share corresponding TCRs, while non-cross-reactive peptides are randomly selected from\npeptides that do not bind these TCRs. This framework ensures that the model learns to identify\npeptides likely to elicit distinct T-cell responses (non-cross-reactive, thus potentially immunogenic)\nversus those likely to cross-react (similar to self, thus less immunogenic).\nManuscript updates:\nUpdated Figures: Fig. S1a\nUpdated Results: Line 143–159\nWe implemented a two-stage training process using TCR-pMHC binding data from IEDB48 and VDJdb49\n(Supplementary Data 1) and T-cell assay data from CEDAR50 (Supplementary Data 2). In the first stage,\nwe constructed a cross-reactive peptide triplet dataset from binding data, where each triplet consists of a\nseed peptide, a cross-reactive peptide (peptide sharing TCR binding with the seed), and a non-cross-\nreactive peptide (peptide not sharing TCR binding). We curated 1,153 MHC-I and 261 MHC-II cross-\nreactive pairs and sampled 10 non-cross-reactive peptides for each pair, yielding 11,530 MHC-I and 2,610\nMHC-II triplets for training. The model was trained to minimize the distance between cross-reactive pairs\n(seed and cross-reactive peptides) and maximize the distance between non-cross-reactive pairs (seed and\nnon-cross-reactive peptides) (Fig. S1a, Methods). In the context of a tumor, mutant peptides more similar\nto wild-type (cross-reactive with the same TCRs) are less likely to be immunogenic because T cells\nrecognizing these self-similar epitopes are typically eliminated during central tolerance. In contrast, more\ndissimilar mutants (less cross-reactive) are more likely to be immunogenic, as T cells capable of recognizing\nthese neoantigens escape thymic negative selection and remain available in the peripheral repertoire.\nNotably, the geometric distance calculated by NeoPrecis-Immuno outperformed the BLOSUM62 and\nPMBEC (an amino acid similarity matrix for peptide-MHC binding)51 distances in differentiating cross-\nreactive and non-cross-reactive peptide pairs (Fig. 2b, Fig. S2a).\nUpdated Methods: Line 521–537\nTo train the immunogenicity model to distinguish mutated peptides from wild-type peptides, we constructed\na dataset of cross-reactive peptide triplets. Each triplet consists of a seed peptide (a reference peptide),\na cross-reactive peptide (sharing TCR binding with the seed), and a non-cross-reactive peptide (not sharing\nTCR binding) (Fig. S1a). The model was trained to minimize the distance between seed-cross-reactive\npairs (similar peptides) while maximizing the distance between seed-non-cross-reactive pairs (dissimilar\npeptides). To mimic the effect of single amino acid substitutions, the difference between both cross-reactive\nand non-cross-reactive pairs was constrained to one Hamming distance. Cross-reactive pairs were\nsampled from each peptide-MHC complex group to ensure functional similarity between the seed and\ncross-reactive peptides. All possible pairs were sampled, yielding 1,153 MHC-I and 261 MHC-II cross-\nreactive pairs.\nFor each cross-reactive pair, non-cross-reactive peptides were randomly selected to be one Hamming\ndistance away from the seed but absent from the cross-reactive set, which comprised peptides within the\nsame peptide-MHC complex group. A BLOSUM62 matrix was used as a prior, ensuring that the substitution\nbetween seed and non-cross-reactive peptides resulted in a negative BLOSUM62 score, indicating\nsubstantial sequence dissimilarity. To enhance dataset diversity, we applied 10-fold generation, ultimately\nproducing 11,530 MHC-I and 2,610 MHC-II triplets (Supplementary Data 1).\n2. It would be appreciated if the authors labeled the training dataset sizes for the models (e.g., how\nmany seed peptides were used from IEDB/VDJdb, and how many data points in CEDAR were\nretained for training, given that it was also used for internal validation).\nWe thank the reviewer for this request. This information was previously provided in the Methods\nsection, and we have now added key sample sizes to the Results section as well for improved clarity\nand accessibility.\nCross-reactive peptide triplet dataset (pre-training):\nThe dataset comprises 1,153 MHC-I and 261 MHC-II seed peptides. For each seed peptide, we\nidentified one cross-reactive peptide (sharing the same TCR) and sampled 10 non-cross-reactive\npeptides. In total, the pre-training dataset contains 11,530 MHC-I triplets and 2,610 MHC-II triplets.\nCEDAR immunogenicity dataset (fine-tuning):\nAfter filtering for predicted MHC binding (netMHCpan %rank ≤ 2 for MHC-I and ≤ 10 for MHC-II), the\ndataset comprised 4,176 MHC-I and 87 MHC-II samples. These were split into training (75%),\nvalidation (10%), and testing (15%) sets, yielding 3,143 training, 415 validation, and 618 testing\nsamples for MHC-I, and 54 training, 11 validation, and 22 testing samples for MHC-II. From the initial\ntesting set of 640 samples, we retained 438 unique peptides that were confirmed to be absent from\nthe training data of all comparison methods (NeoPrecis-Immuno, PRIME, DeepNeo, and ICERFIRE)\nfor unbiased evaluation.\nManuscript updates:\nUpdated Results:\nLine 147–149:\nWe curated 1,153 MHC-I and 261 MHC-II cross-reactive pairs and sampled 10 non-cross-reactive peptides\nfor each pair, yielding 11,530 MHC-I and 2,610 MHC-II triplets for training.\nLine 160–165:\nIn the second stage, we fine-tuned the model using T-cell assay data with T-cell activation labels, enabling\nit to learn the immunogenicity of target peptides (Methods). We curated 4,176 MHC-I and 87 MHC-II\nsamples predicted to bind MHC (NetMHCpan %rank ≤ 2 for MHC-I and ≤ 10 for MHC-II) from CEDAR and\nsplit them into training (75%), validation (10%), and testing (15%) sets. The resulting immunogenicity model\nwas subsequently validated on both an internal and an external testing set to benchmark its performance.\n3. For evaluation metrics, could the authors also report PPVn in the context of neoantigen\nimmunogenicity prioritization?\nThank you for the recommendation. We add PPV@k for the neoantigen immunogenicity prioritization,\nespecially for the NCI dataset. We select k as the number of true positives.\nManuscript updates:\nUpdated Figures: Fig. S2e (top), Fig. S6b (bottom)\nUpdated Results: Line 182–185\nFor MHC-I, NeoPrecis-Immuno outperformed all other predictors across AUROC, AUPRC, and PPV\nmetrics. For MHC-II, NeoPrecis-Immuno achieved the highest AUROC while performing comparably to\nDeepNeo in AUPRC and PPV (Fig. 2d, Fig. S2e).\n4. How are the core positions for Class II peptides determined?\nCore binding regions for both MHC-I and MHC-II peptides are determined using NetMHCpan\npredictions. In addition to binding affinity and rank scores, NetMHCpan outputs the predicted 9-mer\ncore binding region for each peptide-MHC interaction. This 9-mer core is used as input to our model\nfor both MHC class I and class II peptides, ensuring consistent representation across alleles and\npeptide lengths. For peptides shorter than 9 residues, the sequence is padded to obtain the 9-mer\nbinding core as specified by NetMHCpan.\n5. In Figure 3D, could the authors comment on the outlier and why motif-embedding did not improve\nperformance in this case?\nWe thank the reviewer for highlighting this outlier, which corresponds to allele A*68:02 at mutant\nposition 3. This case involves only one positive and one negative sample, detailed in the table below.\nLabel WT_aa MT_aa SubDist GeoDist CRD %rank immunogenicity\nNegative S I 2.90 2.79 2.79 1.2 0.39\nPositive N I 5.56 1.47 1.47 0.13 0.48\nIn the substitution distance (SubDist), the N to I substitution (positive sample) has a larger distance\n(5.56) than the S to I substitution (negative sample, 2.90). However, after motif enrichment, the\ngeometric distance (GeoDist) reverses this ranking. This occurs because the motif enrichment applies\ndifferent scaling factors to embedding dimensions. For A*68:02 position 3, the learned scaling factors\nare 4.55 for dimension 1 and 2.96 for dimension 2 (Fig. 3c). In the embedding space, S and I differ\nprimarily along dimension 1, while N and I differ more along dimension 2 (Fig. 3b). Since dimension\n1 has a larger scaling factor (4.55 vs. 2.96), the S-to-I distance is amplified more than the N-to-I\ndistance, resulting in a larger GeoDist for the negative sample (2.79 vs. 1.47).\nNotably, when MHC binding affinity (%rank) is integrated into the final immunogenicity prediction, the\npositive sample correctly receives a higher score (0.48 vs. 0.39) despite the lower CRD. This case\nillustrates that T-cell cross-reactivity distance is not the sole determinant of immunogenicity—here,\nthe strong MHC binding of the positive sample (%rank = 0.13) compensates for its lower CRD,\nultimately producing the correct prediction. This example demonstrates the importance of integrating\nmultiple signals, as MHC binding can be the dominant factor in certain cases.\n6. In Figure 4C/S5B, there doesn’t appear to be a significant gain when comparing A+R with A+P+R.\nCan the authors comment on why this might be the case? Also, the AUROC values for A+R and\nA+P+R are extremely close — it would be helpful to have the values labeled directly.\nWe thank the reviewer for this observation. We have added metric values directly to the bar plots in\nFig. 4c and Fig. S6b (previously S5b) for easier comparison.\nRegarding the similar performance of A+R and A+P+R, this occurs because NeoPrecis-Immuno\nalready incorporates MHC binding predictions as a covariate during training and inference, effectively\ncapturing presentation information within the recognition metric. Feature importance analysis\nsupports this interpretation: in the A+P+R model, recognition dominates while presentation\ncontributes minimally; in contrast, when recognition is absent (A+P model), presentation becomes an\nimportant predictor (Fig. R10A-B). These results demonstrate that integrating multi-dimensional\nfeatures substantially improves immunogenicity prediction compared to single features alone.\nNeoPrecis-Immuno serves as a comprehensive metric that efficiently captures both T-cell recognition\nand MHC presentation information, explaining why explicit inclusion of presentation features provides\nminimal additional benefit when recognition is already included.\nManuscript updates:\nUpdated Figures: Fig. 4c (figure below), 4d (Fig. R10A), S6b (in minor comment 3), S6c (Fig. R10B)\nUpdated Results: Line 289–302\nTo refine predictions, we integrated the mutation-centric features into a logistic regression model. After\nperforming feature selection, abundance features (A): DNA AF, RNA AF, and RNA expression; presentation\nfeature (P): PHBR; and recognition feature (R): NeoPrecis-Immuno were included in the final integrated\nmulti-dimensional model, NeoPrecis-Integrated (Fig. S6a, Methods). Cross-validation repeated 100 times\non the NCI dataset demonstrated that the A+P+R and A+R achieved similar performance (average AUROC:\n0.848 vs. 0.849 for MHC-I; 0.806 for MHC-II), both substantially outperforming single-feature models (Fig.\n4c, Fig. S6b). This occurs because NeoPrecis-Immuno already incorporates MHC binding predictions as\na covariate, effectively capturing presentation information. Feature importance analysis confirmed that\nrecognition dominates in A+P+R while presentation contributes minimally, whereas presentation becomes\nimportant in A+P models (Fig. 4d, Fig. S6c). These results demonstrate that integrating multi-dimensional\nfeatures improves immunogenicity prediction compared to single features alone, with NeoPrecis-Immuno\nserving as a comprehensive metric that captures both T-cell recognition and presentation information.\nFig. R10. Feature importance in multi-dimensional integrated models. Features were standardized prior to\nmodeling, and importance values were derived from logistic regression coefficients.\n(A) Feature importance in the A+P+R integrated model. (new Fig. 4d)\n(B) Comparison of feature importance across models using different feature combinations, including A+P,\nA+R, and A+P+R. (new Fig. S6c)\n7. In Figure 4D, could the authors also show the feature importance of binding/presentation predictions\nfrom NetMHCpan?\nWe thank the reviewer for this suggestion. NetMHCpan binding predictions are represented in our\nanalysis through PHBR (harmonic mean of MHC binding ranks across alleles), which aggregates\nNetMHCpan predictions to the mutation level. PHBR, as the presentation (P) feature, was already\nincluded in the original Fig. 4d. We have now also added feature importance comparisons between\nA+P and A+P+R models in Fig. R10B (new Fig. S6c). As shown, MHC presentation is an important\npredictor in the A+P model. However, in the A+P+R model, its importance is substantially reduced,\nas the presentation signal is effectively captured within NeoPrecis-Immuno, which incorporates MHC\nbinding predictions as a covariate.\nManuscript updates:\nUpdated Figures: Fig. 4d, Fig. S6c (provided in the previous comment)\nUpdated Results: Line 289–302 (provided in the previous comment)\nReviewer #2 (Remarks on code availability):\nThe code is well documented with examples and usage."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have answered my questions. I have no more questions for the authors.", "reviewer_2": "I appreciate the authors’ thorough responses to my previous comments. I have no further remarks."}, "author_response": ""}]

Decision letter