Ethical statement

This research adheres to all applicable ethical regulations. Patient samples collected in this study were obtained following protocols approved by the institutional regulatory committee (Lausanne University Hospital, CHUV). Written, informed consent was obtained from all patients. In vivo experiments were performed in accordance with Swiss ethical guidelines and under approved licenses (see the next section), ensuring compliance with the 3R (replacement, reduction, refinement) guidelines.

Cancer patient data collection

scRNA-seq/scTCR-seq data from tumors were obtained from internal patients with locally advanced (stage III) or metastatic (stage IV) cutaneous melanoma who had progressed on at least one standard first-line therapy. Tumor samples were obtained following surgery and processed as previously described3,19 for single-cell analysis. Briefly, scRNA-seq and scTCR-seq data were aligned to the GRCh38 reference genome using cellranger count (10X Genomics, v.3.0.1) and vdj (10X Genomics, v.3.1.0), respectively. We also applied TRTpred on external data from ref. 14 (melanoma, n = 4), ref. 11 (GI, n = 5), ref. 8 (n = 1 melanoma, n = 2 breast and n = 12 GI) and ref. 12 (lung, n = 4) collected from the Gene Expression Omnibus using SRA Toolkit (v.3.1.0). All data were processed and annotated following the authors’ guidelines. For the ref. 8 and ref. 12 datasets, neoantigen-specific clones were considered as tumor-reactive while the remaining clones were classified as non-tumor-reactive, adhering to the authors’ ‘closed-world’ assumption8,12. All patients are referenced in Supplementary Table 1.

TCR cloning and tumor reactivity validation

TCRs from the ten internal melanoma patients were annotated following the rationale described in ref. 3. Briefly, TCR antitumor reactivity was interrogated by transferring RNA coding for TCRαβ pairs into recipient activated T cells and Jurkat cell line (TCR/CD3 Jurkat-luc cells (NFAT), Promega, stably transduced with human CD8α/β and TCRα/β CRISPR knockout). Electroporated cells were cocultured with interferon-γ (IFNγ)-treated autologous tumor cells and tumor reactivity was assessed through CD137 upregulation or bioluminescence assay for T cells and Jurkat cells, respectively. From the internal dataset (n = 10 patients), 102 tumor-reactive and 123 non-tumor-reactive CD8+ TCRs were used to build TRTpred and 46 tumor-reactive and 44 non-tumor-reactive CD8+ TCRs were used for the model benchmarking (n = 4 additional patients; Supplementary Table 1). Several TCRs were previously described3. For the dose–response curve, autologous activated T cells electroporated with TCR 1, 3 or 5 were cocultured with tumor cells in IFNγ assay, using precoated 96-well ELISpot plates (Mabtech) as described16.

Statistical models to predict tumor reactivity

Two different approaches were used to predict cell-wise tumor specificity: the signature score and the LR approach. The models first predict cell-wise tumor specificity from scRNA-seq data which is then inferred on the TCR repertoire. The clone-wise score corresponds to the maximum tumor-reactive score obtained by any cell from a given clone.

The signature score approach, standard in RNA sequencing data analyses, uses differential expression analysis methods to derive a signature of tumor specificity which is then used to score cells. To allow comparison between the training and testing scores, the score is scaled based on the mean and standard deviation obtained from the training data. A threshold is then identified to stratify cells into tumor-specific and nonspecific cells, maximizing accuracy. Note that the score thresholds are identified in the training data and applied on the testing data. For this approach, nine different models were constructed upon the selection of nine differential expression analysis methods (DESeq2-Wald/-LRT, edgeR-QFL/-LRT, limma-trend/-voom, Seurat-LR/-negbinomial/-wilcox). These methods were carefully selected based on differential expression analysis benchmarking articles25,26 and applied following best practices (see the section ‘Differential expression analysis methods’). Other parameters such as how to select the genes from the differential expression analysis method (that is, by considering only the log-fold-change or the P value), the number of genes in the signature (that is, signature length), whether to take only the upregulated genes or both up- and downregulated genes (that is, signature side) and the signature score method (Average, AUCell, Singscore, UCell) were defined as hyperparameters to fine-tune (Extended Data Fig. 2a).

The LR approach uses a standard LR coupled with an elastic-net regularization. The type and strength of regularization are defined, respectively, by the α and λ hyperparameters, where α = 0 corresponds to Ridge regression and α = 1 to Lasso regression. From this definition, we derived LR models based on two different feature spaces: first the scaled expression data (RNA) and second the principal components. For the RNA feature space, the genes correlating more than 80% were removed. Because of the large dimensionality, and in addition to the elastic-net, we also tested filtering nonessential features. For both feature spaces, two different dimensionality reduction methodologies were used. The first consists of keeping only features associated with tumor specificity (Wilcoxon, P < 1%). For the RNA feature space, we also applied the differential expression analysis method mentioned above to keep only significant genes (Bonferroni-adjusted P < 0.05). Finally, another model was constructed solely on principal components, explaining more than 10−4 of the explained variances. The combination led to ten RNA- and two principal-components-based models (Extended Data Fig. 2a).

Training and evaluation framework

The evaluations of the 21 models and their associated hyperparameter combinations were performed using an NCV (Extended Data Fig. 2b). This robust framework allows us to iteratively train and test the models on different data partitions called folds. For this application we used an LOPO NCV designed to partition the data into training folds composed of data from all patients but one, which constitutes the testing fold. Iteratively, this approach allows us to simulate the evaluation of the model on new unobserved data from other patients. For the sake of robustness, the performance metric chosen to evaluate the model is the reliable Matthews correlation coefficient (MCC)27. The LOPO evaluation serves as the final model evaluation, and the best model’s hyperparameters are fined tuned using a similar LOPO cross-validation. Ultimately, a final tumor specificity model is obtained by training the model using its best configuration on the whole input data.

Signature collection and comparison

Five CD8+ TIL tumor-reactive signatures from different tumor indications were collected (refs. 8,11,12,14,17) and compared with TRTpred. All signatures have defined an upregulated gene-set but only refs. 14 and 12 have a downregulated gene-set. To compare two signatures, we have used Venn diagrams and the Jaccard index computed as the number of genes in common (intersection) divided by the total number of genes (union).

Model benchmarking

To establish the robustness of the model’s association with tumor reactivity, we subjected the training and evaluation framework to y-randomization tests by applying the same methods with data composed of randomly permuted tumor-reactive annotated clones, repeating the process 100 times. This extensive analysis yielded an average MCC of 0 and an average accuracy of 50%, indicative of the model’s immunity to spurious learning and providing strong validation for its credibility. Encouragingly, we further validated the efficacy of our model on 100 and 205 CD8+ tumor-reactive annotated clonotypes, sourced from four internal and four external (ref. 14) melanoma tumor biopsies. To further explore the generalization potential of TRTpred, we also applied it on three external CD8+ TIL datasets from different tumor indications (ref. 11 (GI, n = 5), ref. 8 (melanoma, n = 1; breast, n = 2; GI, n = 12) and ref. 12 (lung, n = 4)). For the sake of completeness, we collected CD8+ TIL tumor-reactive signatures from these studies and applied them to each dataset. The external signatures were applied on each dataset by using the signature score method described in the respective studies. If not mentioned otherwise, a simple average signature score was computed. Finally, the discriminant power of the different signatures and TRTpred on each dataset was obtained by plotting receiver operating characteristic (ROC) curves and by computing the AUC for the ROC curves.

Differential expression analysis methods

We applied nine different differential expression analysis methods found to work best in scRNA-seq data and applied them following best practices guidelines26,28. These methods were grouped into two classes: the methods developed for scRNA-seq data (LR, negative-binomial and Wilcoxon) and the methods originally developed for bulk RNA sequencing data (edgeR, DESeq2 and limma), named pseudo-bulk methods for the sake of clarity. The single-cell methods were implemented using the Seurat FindMarkers function on the log10-normalized unique molecular identifier counts with all filters (min.pct, only.pos, logfc.threshold, min.cells.group) disabled. To ensure the performance of the pseudo-bulk methods, we applied them on the clone-average log10-normalized unique molecular identifier counts. We did this to obtain a dataset resembling the bulk RNA sequencing data distribution (which reduces inconsistencies in pseudo-bulk methods28) while retaining the behavior of the clones transcription. EdgeR was applied both with the likelihood ratio test (edgeR-LRT, with default dispersion estimate) and with the quasi-likelihood approach (edgeR-QFL). DESeq2 was applied both with the Wald test of the negative-binomial model coefficients (DESeq2-Wald) and with a likelihood ratio test compared with a reduced model (DESeq2-LRT). Limma was applied using two approaches: one incorporating the mean-variance trend into the empirical Bayes procedure at the gene level (limma-trend) and the other incorporating the mean-variance trend by assigning a weight to each individual observation (limma-voom). The log-transformed counts per million values computed by edgeR were provided as input to limma-trend.

Tumor microdissection and RNA extraction

Tumor microdissection and RNA extraction were performed as previously described29. Consecutive sections from fresh-frozen tissue blocks were cut in a cryostat at 8-μm thickness, mounted on precooled PET slides (Leica) at −20 °C for 1 h and fixed in ethanol. They were stained with cresyl violet, cleared in ethanol and microdissected within 20 min after staining using the Leica LMD7000. Laser parameters were set as follows: laser power of 39 mW, a wavelength of 349 nm, pulse frequency of 664 Hz, pulse energy of 58 μJ. Microdissected tissues were collected in 0.5-ml tubes in RNAlater solution (Thermo Fisher Scientific) and kept at −20 °C until RNA extraction. RNA was extracted using RNeasy Plus Micro Kit (Qiagen). To quantify total RNA, Qubit RNA HS Assay Kit and Qubit Fluorometer (Thermo Fisher Scientific) were used. The 2100 Bioanalyzer (Agilent) was used to analyze RNA fragment size using Eukaryote Total RNA Pico assay (Agilent).

Sequential multiplexed immunohistochemistry

Sequential multiplexed immunohistochemistry was performed as previously described29. Fresh-frozen tissue sections were cut at 4 μm, fixed in paraformaldehyde (4%) overnight and permeabilized in 0.5% Triton X-100 in PBS for 30 min. After heat-mediated antigen retrieval in pH 6 citrate buffer for 10 min, endogenous peroxidases, nonspecific proteins, endogenous biotins and avidins were blocked (Dako). After application of the first primary antibody, a biotinylated secondary antibody and a streptavidin-HRP complex were added. Staining was revealed using AEC Chromogen. Tissues were counter-stained using Harris hematoxylin for 1 min and coated by a glass coverslip using an aqueous mounting solution. Slides were scanned into MRXS images using a Pannoramic 250 Flash III scanner (3D Histech). Glass coverslips were removed by immersion in hot water. AEC staining removed by immersion in ethanol of increasing concentrations. Antibodies were stripped by boiling tissue sections in a solution of citrate buffer (pH 6) for 10 to 20 min. Putative residual antibodies were blocked with Fab fragments. Multiplexed immunohistochemistry consisted of sequential cycles of: (1) staining with primary antibodies revealed by AEC Chromogen; (2) tissue section scanning; (3) removal of AEC Chromogen with ethanol; and (4) antibody stripping and blocking with Fab fragments. Primary antibodies were FOXP3 (clone ab99963, Abcam, dilution 1:50) and CD8 (clone C8/144B, Dako, dilution 1:20).

Bulk TCR α and β sequencing

Bulk TCR sequencing analyses were performed as previously described30. Briefly, messenger RNA was isolated and amplified by in vitro transcription. A 5′ adapter was added by multiplex reverse transcription and TCRs were amplified using one primer in the adapter and one in the constant region. Libraries were sequenced on an Illumina instrument and TCR sequences processed using an ad hoc Perl script.

TCR repertoire metrics

To analyze TCR repertoires, two metrics were used: (1) the richness corresponding to the total number of unique clones present in the repertoire and (2) the clonality described by the metric 1-Pielou’s evenness, as previously described31.

MixTRTpred—integration of TCR structural avidity and TCR clustering

Predictions of TCR structural avidity were performed as described32. In brief, a binary LR, based on the CDR3β amino acids that are sufficiently solvent-exposed to interact with the cognate peptide, was used to determine whether a TCR was likely to bind the corresponding pMHC with high or low structural avidity (koff). Avidity levels were also computed with this model on assessable patients, that is, patients with scTCR-seq data composed of both alpha and beta CDR3 chain information. TCR clustering using TCRpcDist was performed as described18. TCRpcDist is a novel and fast structure-based approach that calculates similarities between TCRs using a metric related to the physicochemical properties of solvent-exposed amino acids of the most important residues of this receptor.

The TIL repertoire of patient 14 was analyzed and filtered through the three predictors, TRTpred, the high structural avidity predictor32 and TCRpcDist18. To combine the three axes, we first filtered out the low structural avidity predicted clones and ranked the resulting clones according to their tumor-reactive score. Finally, TCR clustering was applied on the subset of the top 20 tumor-reactive clones. The distance matrix obtained through TCRpcDist18 went through hierarchical clustering (agglomerative method: unweighted pair group method with arithmetic mean), leading to a dendrogram. We chose to split the latter into five distinct TCR clusters given the downstream in vitro and in vivo validation. The selection of five clusters was arbitrary and can be adapted depending on the clinical context.

In vivo study

The in vivo study was performed as previously described16 and was approved by the Veterinary Authority of the Canton de Vaud (under license 3746) and performed in accordance with Swiss ethical guidelines. In brief, Interleukin-2 (IL-2) NOG mice (Taconic Biosciences) were monitored three times a week and given a score based on their weight, behavior, physical condition, dehydration, breathing and tumor burden. As per the protocol, animals reaching a defined score were killed.

TCR transduction for in vivo experiment

TCRs 1, 3 and 5 were selected among the five in vitro-validated tumor-reactive TCRs from patient 14 (Supplementary Table 6) to be tested in vivo. The TCR transduction for the in vivo experiments was performed as previously described16. Briefly, primary CD8+ T cells from a healthy donor were negatively selected with beads (Miltenyi Biotec), activated and transduced as previously reported16,33. Transduced cells were stained with an APC-conjugated anti-mouse constant beta antibody (eBioscience), followed by sorting with anti-APC microbeads (Miltenyi Biotec). Sorted TCR-transduced CD8+ T cells were then expanded and tumor reactivity was assessed in IFNγ ELISpot assays (Mabtech). Transduced cells were then plated at 5 × 103 cells per well and challenged with IFNγ-treated autologous tumor cells at a 1:1 ratio. After 18–20 h of incubation, cells were removed, the plate was developed according to the manufacturer’s instructions and cells were counted using a Bioreader 6000-E (BioSys).

Adoptive T cell transfer in immunodeficient IL-2 NOG mice

The adoptive T cell transfer in immunodeficient IL-2 NOG mice was performed as previously described16. IL-2 NOG mice were subcutaneously injected with 106 autologous human melanoma tumor cells from patient 14 and, once tumors became palpable (day 11), 1–5 × 106 human tumor-specific CD8+ T cell clones were injected in the tail vein, according to the treatment arms described in Extended Data Fig. 5c, with 4–6 mice per condition except in some cases where 3 mice were considered.

Data analyses and computation

Data analyses were performed using R Statistical Software (v.4.0.3). All data processing and analysis was performed using the R dplyr (v.1.1.0) and base libraries. The nested and simple cross-validations were performed using an in-house R library developed to control the models and hyperparameters throughout the folds. The R library glmnet (v.4.1-6) was used to build the LR models and their specifications. Parallelization of the computation was allowed using the foreach (v.1.5.2) library. The differential expression analysis methods were computed using the appropriate R libraries (Seurat (v.4.3.0), limma (v.3.50.3), edgeR (v.3.36.0), DESeq2 (v.1.34.0)) as well as the signature score methods (AUCell (v.1.16.0), UCell (v.1.3.1), singscore (v.1.14.0), GSEABase (v.1.56.0)). Statistical analyses were performed using the standard stats (v.4.1.2) library. The statistical tests used and their specifications are described in the figure legends. Parametric tests, for comparing two or more groups, were applied only on normally distributed variables validated with Anderson–Darling, D’Agostino–Pearson omnibus, Shapiro–Wilk and Kolmogorov–Smirnov tests (GraphPad v.9.1.0); otherwise, nonparametric tests were used.

Plotting description

The figures were generated in R Statistical Software (v.4.0.3) with the ggplot2 (v.3.4.4) R package. Alluvial plots were generated using the ggalluvial (v.0.12.5) R package. The distance heatmaps were performed using the pheatmap function from the pheatmap (v.1.0.12) R package. Plotting of the scRNA-seq-derived UMAP was achieved using Seurat (v.4.3.0) R package functions. Venn diagrams were obtained using the ggven (v.0.1.10) R library. Schematic figures of T cells, cancer cells, TCRs, skin, lungs, intestine, breast, tumor and stromal microenvironment, plate and mouse in Fig. 1a were adapted from templates on BioRender.com. All figures were reprocessed using Adobe Illustrator 2023 (v.27.9.1) solely for esthetical purposes.

Reporting summary

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