DatasetThe limitation of positive TCR–pHLA binding data often motivates data aggregation from multiple sources to build a more substantial dataset for training and evaluation36,55. We collected positive binding triples for training from three publicly available datasets: McPAS-TCR56, VDJdb57 and pMTnet36 (Supplementary Table 1). Through data curation, we retained only those triples specific to Homo sapiens, HLA class I and TCRs that feature only the CDR3 β-chain, as these are critical for determining antigen binding specificity. We also excluded records in VDJdb57 with 0 confidence score. After data preprocessing (Supplementary Note 5), 32,508 unique TCR–antigen–HLA binding triples were obtained for 607 antigens presented by 65 HLA-I molecules and 29,687 TCRs.We acquired two independent external data for testing with strict quality control and standardized preprocessing: (1) 489 experimentally validated TCR–antigen–HLA binding triples from pMTnet36, collected from 25 published works subject to systematic validation by those prior studies (Supplementary Table 1); (2) 425 binding triples from a series of studies of melanoma, lung cancer, head and neck squamous cell cancer, lymphoma and GBM (Supplementary Table 1). These sources detected T cell activation via specific pHLA58,59. Among them, 72% were based on peptide–HLA multimers, 8% were obtained through surface plasmon resonance and the remaining 20% were based on in vitro functional assays (CD137/4-1BB flow cytometry, IFN-γ ELISpot and IFN-γ ELISA). All the CDR3β sequences were acquired through a TCR sequencing assay. The data were subject to strict quality control by unifying the naming conventions and eliminating sequences that are incomplete or contain non-standard amino acids (Supplementary Note 5).To rigorously evaluate the generalization capacity of different models on new sequences, we excluded all the test triples whose antigen–HLA pairs were previously encountered in the training dataset. Additionally, we explored several negative sampling schemes recommended in systematic studies28,30,32,60, including the following. (1) Randomly shuffled sequence triples in the positive data as negative samples31,33,35,36. (2) Unified epitope negative sampling in which the epitopes are sampled by their frequency distributions in the positive dataset32. (3) Reference TCR negative sampling in which each epitope is combined with TCRs sampled uniformly from a reference TCR dataset collected from healthy donors61,62 in which all the TCRs were exposed to all the tested pHLA multimers and no binding signals were detected34,63. In our experiments, we generated negative samples that are ten times larger than the positive ones.PISTEIn the Transformer model40, the combination of positional encoding and semantic embedding leads to attention scores that no longer provide meaningful estimations of token relationships64, particularly when dealing with two or more interacting sequences (Supplementary Note 2).To solve this problem, we propose PISTE. Our intuition is that residues typically move along the cumulative forces acting on them due to their interactions before reaching a stable conformation. Leveraging this insight, we use the attention map in a Transformer as a conceptually appealing alternative to quantify pairwise residue interactions, which then serves as the driving force to update residue positions in an iterative and coherent manner. This strategy combines the flexibility of the Transformer with the consistency of physical priors, leading to useful features for predicting biosequence interaction.The network is shown in Fig. 1 and the three basic building blocks are discussed below.Sequence encoder moduleWe use one-dimensional (1D) convolution to encode the local and shift-invariant features from TCR, peptide and HLA sequences to capture useful and transferable sequence information from short amino acid segments. Three convolutional layers are adopted with a kernel size of 1 × 3, a stride of 1 and skip connections. Using PyTorch’s nn.embedding function, we randomly initialized 64-dimensional vectors to represent 21 amino acid types. These embeddings are updated through backpropagation during training.Sliding-attention moduleSliding attention is a physics-inspired dynamic process that steers the positioning of the residues along the gradient field of their interactions. In this process, the attention (or interaction) between two residues takes into account both their spatial proximity and featural correlations. Then, a series of mode-seeking iterations are used to iteratively ‘drag’ the residues in one sequence towards those of another sequence based on the magnitude of residue interactions (attention). This process allows two or more sequences to virtually ‘slide’ against each other in search of potentially the most stable binding configuration.Sliding attention is defined for two sequences U = {u1, u2…um} and V = {v1, v2…vn}, where ui is the ith residue in U and vj the jth residue in V. We treat V as the reference sequence and U as the sliding sequence. Two concurrent attention views for U and V are computed as follows.Spatial attention. We use an m × n proximity matrix S whose ijth entry signifies the spatial closeness between ui and vj. Here Sij is parameterized by the relative distance between residues ui and vj. For the reference sequence V, its residue positions QV = [q1 q2…qn] are constant integers from 1 to n to signify the linear chain structure of the sequence. For the sliding sequence U, its residue positions are a series of real variables PU = [p1 p2…pm] that are fully optimizable to recover the spatial relations between the residues in U and V. A Gaussian function \({\mathbb{g}}\) is used to estimate Sij as$${S}_{ij}={\mathbb{g}}({p}_{i},{q}_{j})=\exp \left(\frac{-| {p}_{i}-{q}_{j}{| }^{2}}{2{h}^{2}}\right).$$
(1)
Featural attention. We adopt an m × n affinity matrix A whose ijth entry reflects the tendency of two residues ui and vj to interact based on their respective embedding vectors by a function \({\mathbb{f}}\):$${A}_{ij}={\mathbb{f}}({\bf{x}}({u}_{i}),{\bf{x}}({v}_{j}))=\exp \left(\frac{{[{E}_{{\rm{S}}}{\bf{x}}({u}_{i})]}^{\top }[{E}_{{\rm{R}}}{\bf{x}}({v}_{j})]}{\sqrt{d}}\right).$$
(2)
Here x(⋅) is a function that converts a discrete residue type to a d-dimensional vector. The exponentiated inner product is used to estimate the non-negative affinity between two residues, where ES and ER represent the learnable transform matrices for the sliding sequence and reference sequence, respectively.We combine the two attention views and use the non-negative, multiplicative term Wij = Aij ⋅ Sij as a comprehensive indicator of whether residue ui and vj are likely to interact, that is, they have to be both spatially close and exhibit a high affinity to form a strong contact. We further use a 1/0 mask function \({\mathbb{M}}\) to refine Aij values by \({\mathbb{M}}({A}_{ij})={M}_{ij}\cdot {A}_{ij}\) to emphasize the residue pairs forming hydrogen bonds, ionic bonds or hydrophobic interactions (Supplementary Note 11).Using these definitions, we can establish an iterative process to systematically update the positioning of residues (ui) in the sliding sequence U based on their interactions with the reference sequence V. The residue ui location (pi, for i = 1, 2…m) are updated as follows:$$\begin{array}{l}{p}_{i}^{(t+1)}=\frac{\mathop{\sum }\nolimits_{j = 1}^{n}{\mathbb{M}}\left({A}_{ij}\cdot {S}_{ij}^{(t)}\right)\cdot {q}_{j}}{\mathop{\sum }\nolimits_{j = 1}^{n}{\mathbb{M}}\left({A}_{ij}\cdot {S}_{ij}^{(t)}\right)},\\ \,\text{s.t.}\,\,\,\,{S}_{ij}^{(t)}=\exp \left(\frac{-| {p}_{i}^{(t)}-{q}_{j}{| }^{2}}{2{h}^{2}}\right).\end{array}$$
(3)
Here the superscript t is the number of iterations. Considering that Aij is only dependent of the semantic embedding of the residues and is a constant with respect to the residue locations \({p}_{i}^{(t)}\), and that \({S}_{ij}^{(t)}\) is a Gaussian kernel evaluated on the distance between a pair of residues, equation (3) is very similar to the mean shift mode seeking65, for which it has been shown that one such iteration is actually a move (of the residue location \({p}_{i}^{(t)}\)) along the gradient of an underlying density function \({\mathbb{D}}(\cdot )\) with adaptive step size. In our context, this ‘density function’ is the accumulated magnitude of the interactions that the residue ui receives when it is located at position \({p}_{i}^{(t)}\), as \({\mathbb{D}}(\,{p}_{i}^{(t)})=\mathop{\sum }\nolimits_{j = 1}^{n}{\mathbb{M}}\left({A}_{ij}\cdot {S}_{ij}^{(t)}\right)\).It is noteworthy that the positional shift in residue ui due to equation (3) is along the direction of accumulated attractions that residue ui receives at location \({p}_{i}^{(t)}\), by noting \({p}_{i}^{(t+1)}-{p}_{i}^{(t)}=\frac{\mathop{\sum }\nolimits_{j = 1}^{n}{\mathbb{M}}\left({A}_{ij}\cdot {S}_{ij}^{(t)}\right)\cdot \left({q}_{j}-{p}_{i}^{(t)}\right)}{\mathop{\sum }\nolimits_{j = 1}^{n}{\mathbb{M}}\left({A}_{ij}\cdot {S}_{ij}^{(t)}\right)}\); here \({\mathbb{M}}\left({A}_{ij}\cdot {S}_{ij}^{(t)}\right)\) conceptually signifies the magnitude of attraction between residue ui and vj at step t, and \(\left({q}_{j}-{p}_{i}^{(t)}\right)\) signifies the direction of attraction pointing from ui to vj at step t. The bandwidth h in equation (1) controls the size of the receptive field: a larger h allows ui to be attracted to more distant residues in the reference sequence V.As the iteration continues, residue ui moves along the reference sequence V until reaching a local maximum of the interaction density or moves for a pre-defined number of steps (two–five steps). The mode-seeking iteration in equation (3) allows injecting useful physical prior by incrementally adjusting a residue’s position to increase its interaction, or attention, with residues from a counterpart sequence. Compared with learnable positional vectors40 that are merely updated thorough gradient, our positional variables are structurally constrained and physically more interpretable.At the end of the sliding process, the m × n hybrid attention matrix$$W={\mathbb{M}}\left(A\odot S\right),$$
(4)
will serve as a comprehensive estimation of residue-level interactions. Note that W is unnormalized. Depending on whether it is normalized by rows or columns, we can update the representations for both U and V in the form of cross-attention as$${\tilde{{\bf{X}}}}_{U}\leftarrow {D}_{W}^{-1}W{\,{\bf{X}}}_{V}{E}_{V}+{{\bf{X}}}_{U},$$
(5)
$${\tilde{{\bf{X}}}}_{V}\leftarrow {D}_{{W}^{\top }}^{-1}{W}^{\top }{{\bf{X}}}_{U}{E}_{U}+{{\bf{X}}}_{V}.$$
(6)
Here DW and \({D}_{{W}^{\top }}\) are row-wise and column-wise degree matrices for normalization, XU = [x(u1) x(u2)…x(um)] and XV = [x(v1) x(v2)…x(vn)] are residue embedding matrices for U and V, and EV and EU are linear matrices to turn XV and XU into ‘values’, respectively. In equation (5), U is the query and sequence V is the key; in equation (6), V is the query and U is the key. One can also use the additive version of W in equation (4) as \(W={\mathbb{M}}\left(A+S\right)\), which gives a denser attention matrix than the Hadamard product. Finally, no self-attention is used within each sequence before cross-attention.The sliding attention for two sequences is summarized in Algorithm 1 and illustrated in Supplementary Fig. 2b. A comparison with standard cross-attention is shown in Supplementary Fig. 3.
Algorithm 1
Sliding attention for two sequences.
Input: Sliding sequence U: embedding XU, position PU;
Reference sequence V: embedding XV, position QV;
Learnable parameters: \({E}_{S},{E}_{R},{E}_{U},{E}_{V}\in {{\mathbb{R}}}^{d\times d}\);
Hyper-parameters: mask \({\mathbb{M}}\), bandwidth h, steps T.
Output: updated residue embedding \({\tilde{{\bf{X}}}}_{U}\), \({\tilde{{\bf{X}}}}_{V}\).
// initialize variables
1: Initialize XU and XV by random vectors.
2: Set QV as consecutive integers from 1 to n.
3: Initialize PU by m evenly spaced numbers in [1:n].
// compute featural attention A
4: \(A\leftarrow \exp \left({{\bf{X}}}_{U}^{\top }{E}_{S}^{\top }{E}_{R}{{\bf{X}}}_{V}/\sqrt{d}\right)\)—equation (2).
// update spatial attention S, residue position PU
5: for t = 1 to T do
6: \(S\leftarrow \,\text{Gaussian}\,\left({P}_{U},{Q}_{V},h\right)\)—equation (1)
7: \(W\leftarrow {\mathbb{M}}\left(S\odot A\right)\)—equation (4).
8: \({P}_{U}\leftarrow {D}_{W}^{-1}W{Q}_{V}\)—equation (3)
9: end for
// converged attention matrix W
10: S ← Gaussian(PU, QV, h)—equation (1).
11: \(W\leftarrow {\mathbb{M}}\left(S\odot A\right)\)—equation (4).
// update representations of U and V
12: \({\tilde{{\bf{X}}}}_{U}\leftarrow {D}_{W}^{-1}W{{\bf{X}}}_{V}{E}_{V}+{{\bf{X}}}_{U}\)—equation (5).
13: \({\tilde{{\bf{X}}}}_{V}\leftarrow {({D}_{{W}^{\top }})}^{-1}{W}^{\top }{{\bf{X}}}_{U}{E}_{U}+{{\bf{X}}}_{V}\)—equation (6).
14: return \({\tilde{{\bf{X}}}}_{U},{\tilde{{\bf{X}}}}_{V}\).
The attention matrix in equation (4) (or its additive version) can be naturally used to approximate residue-level contact relations between U and V. The nonlinear nature of W allows capturing complex patterns of residue sequences that may curl up in three dimensions, despite the 1D positional variables in sliding attention. We can further augment the sliding attention by extending the 1D positional variables to higher dimensions, enforcing a smoothness constraint to the shift of neighbouring residues, and considering intrasequence residue interactions. These will be studied in our future research.Alignment-based pooling moduleWe propose a systematic way to turn variable-sized biosequences into fixed-length representations, to avoid arbitrary token shift in sequence cutting or padding. Here we exploit a biological prior that HLA sequences have a stable 3D substructure66. In particular, the α-1 and α-2 domains in the α chain of an HLA molecule are connected by a short peptide in the shape of a β sheet, forming a groove that is the key to antigen binding. This allows defining the HLA pseudo-sequence, that is, the part of HLA sequence that is in close contact with the peptide (within 4.0 Å of the peptide), which consists of 34 amino acid residues or positions along the entire HLA molecule23.The HLA pseudo-sequence was used in several studies of HLA–peptide and pHLA–TCR interactions36,37,67. Since the pseudo-sequence has a fixed length, we use it as a skeleton so that both TCRs and antigens can be projected onto it to convert to fixed-length sequences. The alignment is based on the attention matrix in equation (4), which precisely specifies the residue interactions between the two sequences.To project the representation matrix X of a sequence onto a skeleton sequence X0 (HLA pseudo-sequence), we use the attention matrix W between X0 and X as a bridge and left multiply it with X:$${\tilde{{\bf{X}}}}_{\to {{\bf{X}}}_{0}}={D}_{W}^{-1}\cdot W\cdot {\bf{X}}.$$
(7)
Here W is the attention matrix (4) by treating X0 as the reference sequence and X as the sliding sequence, and DW is the row-wise degree matrix of W. The normalized attention matrix \({D}_{W}^{-1}W\) serves as a probabilistic alignment matrix that maps residues from X to those of X0, effectively reshaping X to the same size of the skeleton sequence X0 based on their residue interactions specified by W.Loss function for imbalanced classificationPredicting TCR–antigen–HLA binding requires identifying a small number of truly binding triples from a large repertoire, that is, the positive and negative classes are highly imbalanced. Therefore, we used the following focal loss68:$${\mathscr{L}}(\,{p}_{i})=\left\{\begin{array}{ll}-\alpha {(1-{p}_{i})}^{\gamma }\log (\,{p}_{i})\quad &{\rm{if}}\,\,{y}_{i}=1\\ -(1-\alpha ){p}_{i}^{\gamma }\log (1-{p}_{i})\quad &{\rm{otherwise}}.\end{array}\right.$$
(8)
Here i is the sample index, yi is the class label and pi is the estimated probability for the ith sample to be positive.Sliding transformer for TCR–antigen–HLA binding predictionThe workflow of PISTE for TCR–antigen–HLA binding prediction is as follows (Supplementary Fig. 2a).
1.
Use HLA as the reference, and let peptide slide against it through the sliding-attention module. This allows updating the HLA and peptide representations (pHLA).
2.
Use the pHLA complex as the reference, and let TCR slide along it through the sliding-attention module. This allows simultaneously updating the representations for the TCR and pHLA complex.
3.
Project TCR and peptide representations onto HLA pseudo-sequence by alignment-based pooling.
4.
The representations of TCR, HLA and peptide are passed to a feed-forward layer to make predictions.
These four steps are connected in an end-to-end framework to allow for simultaneous variable optimization. The order of the four steps is biologically meaningful, that is, the peptide–HLA interaction is modelled first before the interaction between the pHLA complex and TCRs. PISTE predicts ternary TCR–antigen–HLA binding, rather than binary (peptide–HLA or peptide–TCR), by using only ternary binding status as the labels. However, if the peptide–HLA binding status was also known, it could be incorporated in training as well.Performance evaluation metricsThe performance was evaluated by AUROC, AUPR and PPVn.In AUROC (TPR versus FPR for a series of threshold values), the true-positive rate (TPR) and false-positive rate (FPR) are computed as$${\rm{TPR}}=\frac{{\rm{TP}}}{{\rm{TP}}+{\rm{FN}}},\,\,\,\,\,\,{\rm{FPR}}=\frac{{\rm{FP}}}{{\rm{TN}}+{\rm{FP}}}.$$Here TP denotes true positive; FN, false negative; TN, true negative; and FP, false positive.In AUPR, the precision and recall are computed by$${\rm{Precision}}=\frac{{\rm{TP}}}{{\rm{TP}}+{\rm{FP}}},\,\,\,\,\,\,{\rm{Recall}}=\frac{{\rm{TP}}}{{\rm{TP}}+{\rm{FN}}}.$$PPVn is the fraction of the top-ranked n prediction triples that are true positives, defined as$${\rm{PPVn}}=\frac{{\rm{TP}}_{n}}{{\rm{TP}}_{n}+{\rm{FP}}_{n}}.$$PPVn is widely used in immunogenicity prediction studies23,24,27. Here n is chosen as the number of true binders in the data, as per ref. 26.Experiment settingsIn training the PISTE, we used the ADAM optimizer with a mini-batch size of 1,024 sequences (triples) and a learning rate of 0.001 with 200 epochs. Each residue type has a dimension d = 64 and is randomly initialized. In the loss function in equation (8), α = 0.75 and γ = 2. Hyper-parameters were chosen as follows. The bandwidth h in equation (1) was fixed as h = 1. The number of iterations t for sliding attention was chosen from {2, 3, 4, 5}, and the best t was determined as the one that leads to the highest evaluation metric (average of AUROC and AUPR) on the validation set, which was chosen as 20% of the training data (the remaining 80% was used for training the model). The codes were written with PyTorch 1.7 and run on a PC with NVIDIA RTX A6000 GPU and 3.70 GHz CPU.Patient specimen collectionThis study was reviewed and approved by the Institutional Review Board of Shanghai Sixth People’s Hospital (declaration 2023-KY-155K). Informed consent was obtained from all the patients and the study strictly adhered to all the institutional ethical regulations. The tumour tissues and peripheral blood samples from eight patients with primary prostate cancer were attained following surgery at the Shanghai Sixth People’s Hospital (see the detailed clinical characteristics listed in Supplementary Table 8). No patients had undergone immunotherapy treatment before surgery. Samples were snap frozen by immediate immersion in liquid nitrogen and stored at –80 °C for next-generation sequencing by the Shanghai Applied Protein Technology. PBMCs were prepared from fresh whole blood by Ficoll–Paque density gradient centrifugation and in 90% foetal bovine serum + 10% dimethyl sulfoxide (DMSO).WES and RNA-seqDNA extraction was executed from both peripheral blood and tumour tissue samples using the QIAamp DNA MiniKit (Qiagen). Quantification of DNA concentrations was carried out using the Qubit 2.0 fluorometer (Invitrogen). The DNA underwent fragmentation into segments measuring 180–280 bp in length, using a Covaris instrument. The preparation of the sequencing libraries and capture of exons were conducted in strict accordance with the manufacturer’s protocol, utilizing the Agilent SureSelect Human All Exon V5/V6 Kit. The captured exons were amplified linearly by polymerase chain reaction and then checked by quantitative polymerase chain reaction. The sequencing procedure was executed on two lanes of the Illumina HiSeq 4000 v. 2 (Pair End 150 bp) platform, strictly adhering to the manufacturer’s guidelines and recommendations set forth by Illumina.The extraction of RNA from fresh tissues was carried out by utilizing a combination of TRIzol reagent and the RNeasy MinElute Cleanup Kit (Invitrogen). The assessment of RNA quality was conducted using a fragment analyser (Agilent Technologies). The TruSeq Stranded Total RNA kit (Illumina) was used for the preparation of sequencing libraries, which were subsequently subjected to 150 bp paired-end sequencing on a HiSeq 4000 sequencer (Illumina).Finally, we obtained whole-exome sequencing (WES) and transcriptome sequencing data of the tumour tissue and exome sequencing data of match normal sample for each patient.Sequencing data processing and immunogenic neoantigen selectionWES information processing. On successful completion of sample sequencing, we leveraged OptiType v. 1.3.5 (ref. 69) to determine the genotypes of patients’ HLA alleles. Meanwhile, we utilized a general mutation calling pipeline to detect somatic variations in the genome70. Trimmomatic v. 0.39 (ref. 71) was used for the WES data quality control. The processed WES data of the tumour and matched blood (as a source of normal germ-line DNA) from each patient were aligned to the reference human genome (hg38) utilizing the Burrows–Wheeler alignment tool v. 0.7.17 (ref. 72). Preprocessing was carried out following the GATK (v. 4.2.0) Best Practices Workflow73 before variant calling. To perform single-nucleotide variant and insertion/deletion mutation calls, MuTect2 (GATK v. 4.2.0)74, VarScan v. 2.3 (ref. 75) and Strelka2 v. 2.9.2 (ref. 76) were utilized. To eliminate false-positive mutations, all the mutations detected with allelic fractions of less than 0.05 or coverage of less than 10× were excluded. Then, all the mutations were annotated by leveraging Ensembl Variant Effect Predictor77. The Quantitative Biomedical Research Center (QBRC) neoantigen calling pipeline was subsequently used to retrieve HLA-I-binding neoantigens of 8–11-mer length from the mutation data70, and the corresponding wild-type sequences of neoantigens were also recorded. A median of 2,617 HLA–peptide complexes per sample were used to combine with TCRs into TCR–antigen–HLA triples and then run through PISTE.RNA information processing. RNA sequencing (RNA-seq) data were aligned to the reference transcriptome (hg38) using Kallisto v. 0.46.0 (ref. 78) to determine the abundance of gene expression levels, quantified as transcripts per kilobase million (TPM).TCR repertoire data. To enhance the diversity and coverage of the patient TCR repertoire, the TCR data for each patient were sourced from two distinct origins. One portion comprised TCR sequences acquired through the analysis of the patient’s WES and RNA-seq data using the MiXCR v. 3.0.13 algorithm79. The other segment was drawn from publicly accessible TCR sequencing data documented in the literature pertaining to prostate cancer patients80.Input to PISTE. The pHLA (with mutant neoantigen or wild-type antigen) and TCR sequences obtained from each patient were then combined to generate all the possible TCR–antigen–HLA triples. These triples were fed into the PISTE model for predicting the binding status for each triple.Peptides ranking. Meticulous screening procedures were taken to select candidates from among thousands of neoantigens. First, we categorized all the predicted binding TCR–antigen–HLA triples from our model by antigen. This categorization allowed us to assess the potential immunogenicity level of each antigen by counting the number of TCRs binding to it. Here we focused on mutated neoantigens that bind with at least 100 TCRs and whose wild-type counterparts do not bind with any TCR. Additionally, considering that a single gene mutation could produce multiple antigens, we excluded those genes (and thus all the mutated neoantigens they produced) whose expression levels were under 5 TPM (ref. 81). After refining our candidate list, we ranked the genes by their expression levels and assessed each gene sequentially from the top of this list. For each gene, we selected the neoantigen with the highest number of binding TCRs as the ‘optimal peptide’. We continued this process until we had selected three to four optimal peptides for a patient, typically requiring probing of two to four highly expressed genes per patient.Peptide synthesisLyophilized peptides for neoantigens were manufactured at ≥95% purity from GenScript. The peptides were verified by high-performance liquid chromatography and stored at –80 °C for testing the T cell reactivity.Expansion of T cells specific to neoantigensPBMCs obtained from patients were used to assess the T cell response to candidate neoantigens in an ex vivo setting. For in vitro pre-stimulation of antigen-specific T cells, PBMCs were thawed and cultured in RPMI 1640 medium (Thermo Fisher, cat. no. A1049101-01) supplemented with 10% foetal bovine serum and 1% penicillin–streptomycin (Thermo Fisher). The cells were stimulated in 96-well cell culture plates at 1.5 × 105 cells per well pulsed with individual neoantigen (2.5 μg ml–1) in the presence of interleukin-2 (20 U ml–1; T&L Biotechnology). Interleukin-2 and peptide were added on days 3, 6 and 8, with the same concentration as before. Here phytohemagglutinin (PHA; 10 μg ml–1) was used as the positive control and DMSO and unrelated peptide as negative controls. Cells were harvested after 10 days post-stimulation; quantification of peptide-specific T cell immune response intensity was conducted with the IFN-γ ELISA and flow cytometry assay.T cell response analysis by IFN-γ ELISA assayIFN-γ secretion of T cells was measured by ELISA using human IFN-γ ELISA kit (Multi Sciences, cat. no. EK180-96). Briefly, 5 × 104 pre-stimulated PBMCs in RPMI 1640 containing 10% foetal bovine serum and 1% penicillin–streptomycin were added to each well of a 96-well plate with a total volume of 150 μl. The cells were subjected to re-stimulation using a peptide concentration of 2.5 μg ml–1 at 37 °C with 5% CO2 for a duration of 24 h; PHA (10 μg ml–1) was used as a positive control and DMSO and unrelated peptide as the negative controls. The concentration of IFN-γ secretion was measured with the EnVision plate reader (PerkinElmer). A positive response was determined when the secretion of IFN-γ greater than 15.63 pg ml–1 and greater than twice the negative control (DMSO and unrelated peptide), according to standard criteria8,82.T cell response analysis by flow cytometryFor ex vivo intracellular cytokine detection, PBMCs were re-stimulated with 5 μg ml–1 peptide or 50 ng ml–1 PMA (YEASEN) and 1 μg ml–1 ionomycin (YEASEN) in complete RPMI 1640 (Thermo Fisher, cat. no. A1049101-01) with 10 μg ml–1 brefeldin A (MKBio) at 37 °C overnight. Subsequently, cells were harvested and resuspended in phosphate-buffered saline (Gibco). After treatment, cells were stained for 30 min at room temperature with a Zombie Aqua Fixable Viability kit (BioLegend, cat. no. 423101), anti-CD3 (clone HIT3a, PerCP, BioLegend, cat. no. 300325) and anti-CD8 (clone SK1, APC, BioLegend, cat. no. 344721). After washing, cells were fixed and permeabilized (Foxp3/Transcription Factor Staining Buffer Set, Thermo Fisher, cat. no. 00-5523-00). Intracellular cytokines were stained with anti-IFN-γ (clone 4S.B3, PE, BioLegend, cat. no. 502508) for 30 min at room temperature. Cells were washed with the fluorescence-activated cell sorting buffer and collected using a ACEA NovoCyte flow cytometer.To assess the formation of specific CD8+ T cells following antigen peptide stimulation, we conducted activation induction markers experiment. Cells were collected after 10 days of pre-stimulation, and then re-stimulated with antigen peptides. Following overnight incubation, cells were harvested and stained with the Zombie Aqua Fixable Viability kit at room temperature for 30 min. Subsequently, the cells were stained with anti-CD3 (PerCP, clone HIT3a, PerCP, BioLegend, cat. no. 300325), anti-CD8 (clone SK1, FITC, BioLegend, cat. no. 344703), anti-CD137 (clone 4B4-1, APC, BioLegend, cat. no. 309809) and anti-CD69 (clone FN50, PE, BioLegend, cat. no. 310905) antibodies. The dilution of all the antibodies was 1:100. After incubation at room temperature for 30 min, cells were resuspended in the fluorescence-activated cell sorting buffer and analysed using ACEA NovoCyte flow cytometer.The gating strategy is shown in Supplementary Fig. 16a. Density maps were drawn for each cell group using ACEA NovoExpress v. 1.6. CD8+ T cell activation is identified when the proliferation percentage of the IFN-γ+ CD8+ population (or the double-positive rate of CD69 and CD137) following antigen peptide stimulation is 20% higher than that of the control group (DMSO). Meanwhile, we also quantified the proportion of CD8+ T cells in CD3+ T cells after antigen peptide stimulation.Graphical and statistical analysesPlots and analyses were generated using matplotlib and seaborn package in Python v. 3.8; survival package and survminer package in R v. 4.2.2; and GraphPad Prism software v. 8. A two-sided t-test was used to compare the continuous variables between two groups. To accommodate multiple comparisons, a standard one-way analysis of variance with Dunnett’s test was used. To investigate the existence of a positive ordinal association between the INAL and efficacy of immune therapy, we utilized the Jonckheere–Terpstra test. Survival curves were generated through the Kaplan–Meier method, whereas the log-rank test was used to assess the presence of significant differences between two survival curves. P values of <0.05 were considered to be statistically significant.Reporting summaryFurther information on research design is available in the Nature Portfolio Reporting Summary linked to this article.