# An integrative multimodal machine learning signature of primary resistance to immunotherapy in advanced non-small cell lung cancer: biomarker analysis from the PIONeeR study

Fabrice Barlesi, Florence Monville, Laurent Greillier, Natalie Ngoi, Joseph

Ciccolini, Stephane Garcia, Jean-Philippe Dales, Florence Sabatier, L.

Arnaud, Amélie Pouchin, et al.

► To cite this version:

Fabrice Barlesi, Florence Monville, Laurent Greillier, Natalie Ngoi, Joseph Ciccolini, et al. An integrative multimodal machine learning signature of primary resistance to immunotherapy in advanced non-small cell lung cancer: biomarker analysis from the PIONeeR study. 2025. hal-05241459

HAL Id: hal-05241459

https://inria.hal.science/hal-05241459v1

Preprint submitted on 16 Jan 2026

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# An integrative multimodal machine learning signature of primary resistance to immunotherapy in advanced non-small cell lung cancer: biomarker analysis from the PIONeeR study

Fabrice Barlesi, Florence Monville, Laurent Greillier, Natalie Ngoi, Joseph Ciccolini, Florence Sabatier, Stephane Garcia, Jean-Philippe Dales, Cyril Foa, Laurent Arnaud, Amélie Pouchin, Frédéric Vely, Sivan Bokobza, Anastasiia Bakhmach, Andrea Vaglio, Mélanie Karlsen, Paul Dufossé, Clarisse Audigier-Valette, Jacques Letreut, Lilian Laborde, Pierre Milpied, David Pérol, Mohamed Boussena, Celestin Bigarre, Mourad Hamimed, Richard Malkoun, Vanina Leca, Marcellin Landri, Maryannick Le Ray, Marie Roumieux, Julien Mazières, Maurice Pérol, Jacques Fieschi-Meric, Eric Vivier, Sébastien Benzekry

## Highlights:

* We conducted a prospective, multicenter biomarker study of 439 patients with advanced NSCLC treated with anti-PD-(L)1 therapies across first- and later-line settings, profiling a total of 443 clinical, routine laboratory parameters, tumor and blood-derived biomarkers prior to treatment initiation.

* A machine learning-derived 18-feature multimodal signature predicted primary resistance (PrR) to anti-PD-(L)1 therapy with higher accuracy than standard biomarkers (AUC 0.73, PPV 60%). Routine blood test features, including serum chloride, albumin, C-reactive protein, and monocyte-to-lymphocyte ratio, emerged as dominant contributors to the predictive model.

* The model demonstrated robust prognostic stratification for PrR and progression-free survival.

* Standard-of-care immunotherapy biomarkers PD-L1 expression and TMB showed limited predictive utility (AUC ~0.60) and were consistently outperformed by the integrative multivariable model.

* An interactive digital tool using SHAP values enabled individualized prediction interpretation and illustrated heterogeneous, patient-specific resistance drivers.
# ABSTRACT

## Background

Immune checkpoint inhibitors (ICIs) have transformed the treatment landscape for advanced non-small cell lung cancer (NSCLC), yet primary resistance remains common, with only ~50% of patients responding to first-line chemo-immunotherapy and 20-30% to monotherapy. Existing biomarkers such as PD-L1 expression and Tumor Mutational Burden (TMB) demonstrate limited predictive accuracy, underscoring the need for more comprehensive, integrative approaches.

## Methods

We conducted a prospective, multicenter study involving 439 patients with advanced NSCLC treated with anti-PD-(L)1 ICI across first-line combo with chemotherapy and later-line monotherapy settings. A total of 443 pre-treatment tumor and blood-derived biomarkers—including genomic alterations, immune cell phenotypes, proteic markers, and routine laboratory tests—were profiled. Extensive biostatistics adjusted for PD-L1 expression were conducted. A rigorously benchmarked machine learning (ML) pipeline including 36 feature selection methods embedded into an optimism-correction framework was applied to identify predictors of primary resistance (PrR).

## Results

Single biomarkers showed limited predictive utility, with PD-L1 (AUC 0.62, positive predictive value (PPV) 49.6%), TMB (AUC 0.55, PPV 43.1%), and key gene mutations (e.g., STK11, KEAP1) failing to achieve significance after multiple testing correction. A gradient boosting ML model integrating 18 selected features yielded a corrected AUC of 0.69 and a Positive Predictive Value (PPV) of 60% for PrR, outperforming standard biomarkers. In first-line patients, the model achieved a PPV of 51% and Negative Predictive Value (NPV) of 79% (baseline PrR rate: 29.9%); in subsequent-line patients, PPV reached 64% (PrR rate: 55.1%). Importantly, the signature also stratified Progression-Free Survival (PFS): high-risk patients had a median PFS of 3.9 vs. 14.6 months in low-risk patients (HR 0.307, p < 0.0001). Features from routine blood tests— such as serum chloride, albumin, CRP, and monocyte-to-lymphocyte ratio (MLR)— accounted for half of the final model and demonstrated independent associations with both PrR and PFS (e.g., chloride: OR 0.616, AUC 0.626; HR 0.685, C-index 0.61).
SHAP-based individual-level model explainability revealed heterogeneous and nonlinear biomarker contributions, including cases where high CRP, low albumin, or elevated MLR overrode favorable PD-L1 or Treg profiles. A biomarker dashboard including interactive visualizations is available at https://compo.inria.fr/pioneer-website/.

## Conclusions

Multimodal machine learning integration of clinical, genomic, immune, and laboratory data enables improved prediction of ICI resistance in NSCLC beyond current biomarkers. This approach not only captures the multifaceted nature of tumour-host interactions but also highlights the underrecognized predictive value of accessible blood-based markers, offering a path toward individualized immunotherapy decision-making.

**Keywords:** multimodal biomarkers; machine learning; immunotherapy resistance; biomarker discovery
INTRODUCTION

Immune checkpoint inhibitors (ICI) directed against the anti-PD(L)1 axis have led to a paradigm shift in the treatment of patients with advanced wild-type non-small cell lung cancer (NSCLC)¹,². Despite these advances, just 50% and 20-30% of patients respond to combination chemo-immunotherapy or anti-PD(L)1 monotherapy in the first- and second-line treatment setting, respectively³⁻⁹. Identifying nonresponding patients is of primary importance to develop new treatment strategies. While tumoral PD-L1 expression⁶,¹⁰,¹¹ and tumour mutational burden (TMB) ¹²⁻¹⁶ have been widely investigated and approved as biomarkers to predict response to ICI in advanced NSCLC, their clinical utility is constrained by assay variability, limited sensitivity and specificity, and spatiotemporal tumour heterogeneity². Beyond conventional markers, studies have proposed other potential response biomarkers¹⁷,¹⁸ such as HLA heterogeneity¹⁹, T cell receptor repertoire²⁰,²¹, tumour-infiltrating lymphocytes²²,²³, inflammatory¹⁴ or interferon-γ signatures²⁴,²⁵ or the presence of tertiary lymphoid structures²⁶. Yet, limited overlap exists between these individual markers¹³, further underscoring their independent nature and the limitations of single biomarker assessment.

While the majority of biomarker discovery efforts have focused on predicting ICI response, the identification of reliable biomarkers for primary resistance (PrR) remains underdeveloped. Mechanistically, studies have implicated multiple pathways such as deficient antigen presentation (e.g., loss of β₂-microglobulin)²⁷, genomic mutations (e.g. in JAK1/2²⁸, KEAP1 or STK11²⁹⁻³¹) or chromosomal loss (e.g. 9p21 deletion³²) leading to defective interferon-γ signalling in mediating anti-PD(L)1 resistance. Additional molecular drivers of resistance such as elevated TROP2 expression³³, or the expression of other immune-inhibitory molecules have been hypothesized³⁴ but lack consistent validation across studies.

Given the complexity of tumour-host interactions, novel aggregate models that integrate multi-parametric data may better capture the complexity of the tumour immune landscape and tumour-host interaction at an individual patient level, and in turn, better predict ICI response and resistance²³,³⁵⁻³⁸. The core premise of this approach is that integrating orthogonal biological inputs—spanning genomics, proteomics, clinical features, and immune phenotyping—has the potential to uncover interactions not captured by unimodal biomarkers. Machine learning approaches are
especially well-suited to extract predictive patterns from such high-dimensional data
and have already demonstrated utility in combining imaging, clinical, and molecular
modalities to improve predictive performance over standard biomarkers alone,
accelerating biomarker discovery³⁵⁻³⁸. A previous study successfully combined routine
baseline clinical data, PD-L1 immunohistochemistry features, TMB, genomic data, and
computed tomography radiomic data, demonstrating an improved predictive ability of a
machine learning-derived multimodal model for anti-PD(L)1 ICI response in NSCLC
compared with unimodal features or standard-of-care approved biomarkers alone³⁶.
While the study only utilized retrospectively obtained data from routine diagnostic
procedures at a single pre-treatment timepoint, it provided proof-of-concept that
existing available data from multiple modalities can be annotated, abstracted and
combined using machine learning methods for next-generation biomarker
development³⁶.

We hypothesized that tumour-host immune interactions are central to anti-PD(L)1 ICI
resistance and sought to conduct biomarker discovery through a rigorous prospective
multi-parametric profiling of 443 distinct standard clinical parameters, tumour and
blood-derived biomarkers from 439 advanced NSCLC patients prior to the receipt of
standard-of-care anti-PD(L)1 ICI, either alone or combined with chemotherapy, in the
first- and subsequent-line settings. Using both traditional biostatistical and machine
learning methods, we constructed a multimodal signature predictive to primary ICI
resistance that had improved predictive ability compared to individual features or
standard-of-care approved biomarkers.
# MATERIALS AND METHODS

The study was conducted in accordance with the Helsinki declaration, French laws and regulations and the International Conference on Harmonization (ICH) E6 Guideline for Good Clinical Practice. The study was approved by the French ethics committee (Comité de Protection des Personnes Ouest II Angers, no. 2018/08) and the French drug and device regulation agency (Agence Nationale de Sécurité du Médicament, no. 2018020500208). Informed consent was obtained from each participant before any study procedure. The study is registered at [ClinicalTrials.gov (NCT03493581)](https://ClinicalTrials.gov).

## Study design

This is a prospective study conducted in 17 centres in France. Key eligibility criteria included patients with histologically confirmed diagnosis of advanced or recurrent NSCLC planned for treatment with standard-of-care anti-PD(L)1 ICI. Eligible patients receiving front/first-line therapy had to be planned for platinum-based chemotherapy and anti-PD(L)1 ICI combination therapy. Eligible patients receiving second or third line therapy had to have progressed on a prior line of platinum-based chemotherapy and be planned for an approved anti-PD-(L)1 ICI as monotherapy. All patients had to have available formalin fixed paraffin embedded (FFPE) archival tumour for analysis, be 18 years of age or older, have Eastern Cooperative Oncology Group performance status (ECOG PS) of 0 or 1, have adequate organ function and be naïve to anti-PD-(L)1 ICI at the time of study entry. Investigator-assessed radiographic tumour assessment was performed every 6 weeks, using Response Evaluation Criteria In Solid Tumours (RECIST) 1.1³⁹. Upon disease progression on ICI, patients entered a 24-week follow-up period and new anti-cancer therapy and survival data were collected (Supplementary Figure 1).

## Endpoints

The primary endpoint was PrR, defined differently between patients receiving monotherapy or combination therapy, according to the 2023 guidelines from the Society for Immunotherapy of Cancer⁴⁰,⁴¹. For monotherapy, PrR was defined as the best overall response (BOR) being progressive disease (PD) or stable disease (SD) lasting more than 6 months from the first dose of ICI (per RECIST 1.1 as assessed by the Investigator). For combination therapy, PrR was defined as progression within 6
months from the first dose of ICI, irrespective of the best response. In other words, a
combination therapy who initially responded but progressed within 6 months was
considered resistant, while the same case in the monotherapy setting was not defined
as PrR⁴⁰,⁴¹. To be eligible for assessment, patients had to have drug exposure for at
least 6 weeks (monotherapy or combination therapy) or progressing disease before 6
weeks (combination therapy only). The secondary endpoint was progression-free
survival (PFS), defined as time from study inclusion to progression of disease or death.

**Clinical parameters and circulating biomarkers**

Clinical parameters were collected two weeks prior to ICI initiation, these included
standard clinical parameters (e.g.: ECOG PS, heart rate, systolic and diastolic blood
pressure) and routine clinical blood tests (e.g., blood counts or biochemistry, see
**Supplementary Table 1** for a full list). Peripheral blood samples were collected at
baseline and at on-treatment timepoints for cell populations and cytokine analysis.
Immune cell profiling was performed using flow cytometry on fresh whole blood to
quantify T cells, B cells, natural killer (NK) cells, and total lymphocytes. Frozen
peripheral blood mononuclear cells (PBMCs) were thawed, stained for viability, and
analysed by flow cytometry. Serum levels of soluble CD25 were quantified via ELISA,
and key cytokines (interferon (IFN)-γ, tumor necrosis factor (TNF), interleukin (IL)-6,
IL-10) were measured using a cytometric bead array. Matched peripheral blood cell
pellets were preserved as germline references for whole-exome sequencing (WES).
Tumour mutational burden (TMB) was calculated based on coding nonsynonymous
variants per megabase of tumour DNA.

Circulating endothelial biomarkers were analysed from platelet-free plasma collected at baseline and at week 12 post-ICI initiation. Circulating endothelial cells (CECs) were isolated using immunomagnetic beads coated with CD146 and enumerated by fluorescence microscopy. Soluble endothelial adhesion molecules (sVCAM-1, sICAM-1, sE-selectin, sP-selectin) and angiogenic factors (angiopoietin-2, VEGF-A, VEGF-R1, thrombospondin-1, soluble tissue factor (sTF), sCD146) were measured using multiplexed Luminex assays and ELISAs. Endothelial extracellular vesicles (EVs) were identified by flow cytometry as Annexin V+/CD31+/CD41-/CD146+ events. Tissue factor (TF) activity of EVs was evaluated using an FXa generation assay adapted from the Chapel Hill model⁴².
# Tissue-based biomarkers

## Immunohistochemistry

Archival tumour specimens were processed into ≤ 22 consecutive 4-µm FFPE sections per block. Immunohistochemistry (IHC) was performed following a prespecified prioritisation schema. Samples containing <100 viable tumour cells on section 1 were considered non-evaluable. PD-L1 expression was assessed using the Immunoscore IC kit clone HDX3 (Veracyte) on either Ventana Benchmark XT (Roche) or an automated Bond RX platform (Leica) Tumour proportion score (TPS) was calculated as the percentage of viable tumour cells with membrane staining. Multiplex IHC was performed using two panels targeting T cell (CD3, CD8, PD-1, LAG-3, TIM-3) and myeloid cell (CD11b, CD14, CD15, HLA-DR, S100A9, LOX-1) markers on Bond RX platform (Leica). Dual IHC staining for CD8-PD-L1 and CD4-FOXP3 (Veracyte SAS) enabled quantification of immune cell densities and spatial proximity within the same section. Sections were stained using Ventana Benchmark XT (Roche) platforms and visualised with OptiView and UltraView detection kits. Slides were digitised at 20× resolution (NanoZoomer XR, Hamamatsu) and analysed using HALO software (Indica Labs) for region annotation, cell segmentation, phenotype classification, and spatial analytics.

## Immunofluorescence

Multiplex immunofluorescence (mIF) was performed using three mIF panels targeting innate immunity, NK/T cells, and B cells and were applied to sequential FFPE sections. Staining was automated using the Bond RX platform (Leica) with the Opal 6-Plex kit (Akoya Biosciences). Slide imaging was performed on the Phenolmager HT platform at 20× magnification using seven spectral filters. Unmixing and image processing were conducted using inForm (Akoya), followed by quantitative analysis in HALO v3.4 (Indica Labs). Tumour and normal-adjacent tissue regions were annotated manually, and the Highplex FL module in HALO v3.4 was used to identify and phenotype individual cells.

## Genomic sequencing

Macro dissection of tumour regions was performed from 1–14 FFPE slides. Coextraction of DNA and RNA from tumour was done using the Quick-DNA/RNA FFPE Miniprep Kit (Zymo Research). Matched germline DNA was extracted from peripheral blood pellets using the DNeasy Blood & Tissue Kit (Qiagen). DNA quality and
concentration were assessed using spectrophotometric (Nanodrop) and fluorometric (Qubit) methods. Whole exome sequencing (WES) libraries were enriched with xGen Exome Hyb Panel (IDT Technologies) and sequenced on an Illumina NovaSeq platform to achieve a mean target coverage of 100× for tumour and 70× for normal tissue. TMB was computed as previously described¹⁶,⁴³.

Total RNA extracted from FFPE sections was quantified and quality-controlled using Qubit RNA HS (Thermo Fisher) and Agilent ScreenTape. Whole transcriptome libraries were prepared using the QuantSeq 3' mRNA-Seq Kit (Lexogen) and sequenced on an Illumina NextSeq 2000 to generate 75-bp single-end reads.

## Data infrastructure to support multimodal integration

Computational pre-processing, classical biostatistics and supervised analyses were performed in a controlled container (docker using the rocker/tidyverse:4.2.2 image) encompassing operating system and *R* and *python* core and package versions (4.2.2 and 3.10.12). Three packages – all open source and publicly available – were internally developed for providing dedicated tools for dedicated analyses: compo.EDA (biostatistics analyses and reporting)⁴⁴, compo.tidyML (automatic machine learning analyses and reporting)⁴⁵ and ROOFS (RObust biOmarker Feature Selectio)⁴⁶. A GitLab continuous integration/continuous deployment pipeline was implemented to ensure robust integrity testing, traceability and reproducibility of the results during interim data integration and code continuous development. It comprised automatic installation of the environment and packages, testing for data integrity and code, processing of data analysis reports and their continuous deployment to (1) the PIONeeR dashboard and (2) the PIONeeR website. These websites were developed with the ultimate goal of enabling data sharing of the data from PIONeeR with the academic community. The PIONeeR dashboard allows users to perform biostatistical comparisons for all individual biomarkers with the specified outcomes across all biomarker types. The PIONeeR website is a central hub for the extensive data accumulated from PIONeeR to enable optimal, centralized sharing of biostatistical analysis results between PIONeeR Consortium partners. A shiny application was also built to visualize a given patient's clinical history and longitudinal follow-up.

## Log-transformation, missing values and collinearity reduction
After visual inspection, multiple variables showed non-gaussian statistical distributions. (e.g., circulating and tissue biomarkers in absolute units of cells / mm³ and cells / mm²). To fit with underlying assumptions of statistical tests (e.g., t-tests), these and others (239/442, **Supplementary Table 1**) were log-transformed when leading to a variable distribution closer to Gaussian [website].

Biomarkers with high rate of missing values (≥ 70%) were removed (106/442, **Supplementary Table 1**).

To limit redundancy among predictors, multicollinearity was assessed using the variance inflation factor (VIF)⁴⁷. Variables with a VIF greater than 5 were iteratively removed, following a stepwise procedure, until all remaining variables exhibited acceptable collinearity levels (VIF ≤ 5). This process yielded a final dataset of 214 features with reduced interdependence.

## Univariable biostatistical analysis

To conduct univariable analysis of biomarker associations with the binary outcome variable of PrR, logistic regression and appropriate statistical tests were employed. Numeric variables were analysed using t-tests (sample size ≥ 30) or Wilcoxon tests (sample size < 30), while categorical variables were analysed using chi-squared (sample size ≥ 50) or Fisher's tests (sample size < 50). Area under the receiver operating characteristic (ROC) curves (AUC) were also computed for each marker.

For analysis of biomarkers and their association with PFS, proportional hazard Cox regression was performed. This was also assessed with the log-rank test on stratified Kaplan-Meier (KM) curves. For each continuous feature, an optimal cut-off was computed using the maximally selected rank statistic, implemented in the *surv cutpoint* function from the *survminer* package. This cut-off was constrained to yield at least 10% of individuals in each group. Harrell's concordance statistics (or C-index) were also computed for all variables.

Importantly, to identify associations independent to PD-L1 expression, for both logistic and Cox regressions were performed adjusting for PD-L1 (i.e., multivariable models including the biomarker and PD-L1 expression as a continuous variable). To consider multiple hypothesis testing due to the large number of biomarkers, adjusted p-values were computed for each type of statistical analysis according to the Benjamini-
Hochberg false discovery rate controlling procedure<sup>48</sup>. These biostatistical analyses were implemented within a dedicated R package named compo.EDA<sup>44</sup>.

## Supervised Machine Learning

### Pre-processing

Multiple imputation methods (e.g., k-nearest neighbours, multiple imputation) were tested but none improved results (i.e., predictive performances) over simple imputation by the median for numeric variables and the mode for categorical ones. Numerical variables were centered and scaled (z-score normalization). These pre-processing steps were encapsulated as a step of the pipelines described below to prevent any potential information leakage between training and testing sets.

### Optimism-correction evaluation framework

The modeling workflows, including pre-processing (imputation, scaling), feature selection, model selection, and hyperparameter tuning, were implemented via ROOFS<sup>46</sup>, an internally developed Python package built on the scikit-learn library<sup>49</sup>. ROOFS relies on the resampling-based .632 optimism correction method for rigorous performance evaluation<sup>50,51</sup>. Briefly, identical modeling pipelines were applied to both the full dataset and each resample. The estimate of a given performance metric $\theta$ was calculated as a weighted sum of the apparent and out-of-bag performance:

$$ \theta_{.632} = (1 - w) \theta_{app} + w \theta_{out}, $$

where $\theta_{app}$ is the performance metric computed on the full dataset, and $\theta_{out}$ is the performance metric computed on external samples (patients not included in the resampled dataset), averaged across resamples. In the standard .632 method, the weight $w=0.632$ corresponds to the approximate number of samples included in a bootstrap. Thus, weighting takes into account the overlap between the bootstrap samples and the original data. Since we used subsampling without replacement to 80% of the original number of patients, we adjusted the weight to 0.8 accordingly.

One hundred subsamples per pipeline were used. The primary evaluation metrics were the .632 corrected AUC and Positive Predictive Value (PPV) for prediction of PrR. Secondary metrics were the .632 corrected balanced accuracy, Negative Predictive Value (NPV), sensitivity, and specificity.

### Feature selection – Signature
Thirty-six feature selection (FS) methods were benchmarked for performance and stability, defined for a given FS method by⁵²:

$$ \text{Stability} = 1 - \frac{\frac{1}{p} \sum_{f=1}^{p} s_f^2}{\frac{k}{p} \left(1 - \frac{k}{p}\right)}, $$

where *p* is the total number of features in the dataset, $s_f^2 = \frac{B}{B-1} \hat{p}_f(1-\hat{p}_f)$ is the unbiased sample variance of the frequency of selection of the feature *f* over the *B* subsamples and $\bar{k}$ is the average number of selected features in a subsampled dataset (thus $\frac{\bar{k}}{B} \left(1 - \frac{\bar{k}}{B}\right)$ is the expected variance of the frequency of selection for a uniformly randomly selected feature).

The embedded methods included Least Absolute Shrinkage and Selection Operator (LASSO⁵³, stability selection⁵⁴, bootstrap LASSO⁵⁵, adaptive Lasso⁵⁶, relaxed Lasso⁵⁷, sparse group Lasso with correlation and Minkowski distances⁵⁸, Stabl LASSO and adaptive Lasso⁵⁹, and RENT⁶⁰. For filter methods, unless otherwise specified, the number of outputted features was fixed at 15. These methods included Fisher score⁶¹, trace ratio⁶², Conditional Infomax Feature Extraction (CIFE)⁶³, Conditional Mutual Information (CMIM)⁶⁴; Double Input Symmetrical Relevance (DISR) filter⁶⁵, Interaction CAPping (ICAP)⁶⁶, Joint Mutual Information (JMI)⁶⁷, Mutual Information Maximization (MIM)⁶⁸, ReliefF⁶⁹, and minimum Redundancy Maximum Relevance (mRMR)⁷⁰, F-score⁶¹, Gini index⁶¹, t-score⁶¹. Hierarchical clustering was utilized by selecting the feature with the highest correlation with the outcome within each cluster. We also employed univariable PFS-Cox or PrR-logistic regression (controlled by PD-L1) or t-test-based filtering with an adjusted p-value threshold of 0.05. Lastly, wrapper methods comprised recursive feature elimination using either logistic regression or random forest⁷¹, recursive forward selection with random forest, and Shapicant⁷².

All the methods were grouped in an internally developed, open-source python package ROOFS, allowing unified benchmarking⁴⁶. The method that exhibited the best performance (AUC) / stability trade-off was filtering on univariable statistics obtained by Cox regression for association with PFS (Benjamini-Hochberg false discovery rate (FDR), $q < 0.01$)⁴⁸.
This feature selection method was also applied for comparison of performances across data modalities, resulting in per-data modality (block) feature subsets.

## Predictive modelling

### Supervised machine learning (ML)

Ten ML models (e.g., logistic regression, random forests, (extreme) gradient boosting, multi-layer perceptrons,...) were evaluated, resulting in 360 variable selection/model pairs. Models were implemented with scikit-learn⁴⁹. Hyperparameter tuning, tested via randomized search, offered negligible performance gains and was omitted in final models to optimize execution time. For the tree-based models in our main analysis, we limited the depth of the trees to 2, increased the minimum number of patients in a leaf node to 40, and set the number of tree estimators to 200. Logistic regression was fitted without regularization. The rest of the models were used with their default parameters.

Exploration and reporting of ML performances using classical training / test and stratified cross-validation splits were performed using an internally developed, open-source, compo.tidyML R package, built on the tidymodels framework⁷³, enabling comprehensive machine learning workflows, performance evaluation, and automated reporting.

### Validation of the machine learning score for PFS stratification

The final pipeline was also assessed for its ability to predict PFS. A stratified train/test split (2/3 – 1/3) was performed by clustering (k-means algorithm, k=5) according to known prognostic clinical variables (PD-L1 expression at the 1% threshold, histology, performance status and sex) and using the cluster as stratification variable. The model fitted on the training set generated a predictive score (probability) of PrR, further employed to stratify PFS. A cut point for this score was determined on the training set as the optimal one for separating Kaplan-Meier PFS curves using the *surv cutpoint* function from the *survminer* R package (minprop = 0.1).

## Explainability

To interpret the importance of each marker in the contribution to prediction of PrR, SHapley Additive exPlanation (SHAP) analysis was conducted, in both the marker and patient dimensions, using the *kernelshap*⁷⁴ and *shapviz* R packages⁷⁵,⁷⁶. For the individual analysis, the model was trained on the training set described above and
predictions and SHAP individual marker importances reported for patients in the test set.

In addition, line-specific decision trees were built on the signature variables, using the *DecisionTreeClassifier* (max_depth = 4, min_sample leaf = 10, ccp alpha=0.01), from the *scikit-learn python* library⁴⁹.
# RESULTS

## Clinical characteristics of included patients who received first line chemo-immunotherapy or later-line anti-PD-(L)1 monotherapy

498 patients were prospectively screened at 17 centers in France between March 2018 to January 2023. 439 patients met inclusion criteria and were included in the biomarker analysis, and 435 patients were evaluable for the primary endpoint of PrR (Table 1). Clinical characteristics of the cohort are presented in Table 1. Of the patients evaluable for PrR, 268 patients received first-line chemo-immunotherapy, and 167 patients received second- or third-line anti-PD(L)1 monotherapy (Table 1). Radiological response was prospectively assessed by investigators by Response Evaluation Criteria In Solid Tumours (RECIST) 1.1 (RECIST 1.1) criteria and the overall response rate to chemotherapy combined with anti-PD(L)1 ICI in the first-line setting was 52.3% [best overall response: 128/269 with partial response (PR), 15/269 with complete response (CR), 76/269 with SD, 9/269 with PD, 41/269 not available]. In the second- or third-line setting, the overall response rate was 20.0% [best overall response: 26/170 with PR, 8/170 with CR, 65/170 with SD, 36/170 with PD, 35/170 not available].

## Biomarker features, collinearity reduction and outcomes

From all the markers collected during the PIONeeR study, we focused here on the pretreatment (baseline) timepoint prior to commencement of anti-PD(L)1 ICI. We prospectively collected clinical variables and standardized measurements from routine laboratory blood tests performed on the date of, or no more than two weeks before, the first dose of anti-PD(L)1 ICI, except for PD-L1 TPS that was analyzed on archived tumor biopsy for line ≥ 2. Additionally, multi-modal deep phenotyping of tumor and blood samples obtained was performed. In total, these procedures yielded a total of 442 individual biomarker features (tumour-based assays: 227 individual biomarkers, blood-based assays: 156 biomarkers, routine laboratory blood tests: 48 biomarkers, clinical parameters: 11 markers) assessed for each included patient at the pretreatment (baseline) timepoint (Figure 1A-B). 394 of these biomarkers were quantified on a continuous scale, while 48 biomarkers were categorical variables. The global rate of missing values was 39%. 106 biomarkers 70% of missing values were removed (Supplementary Table 1). Most of them (101/106) were from multiplex immunofluorescence, leaving 12 of these to the subsequent analysis. The remaining 5
were the Shannon entropy of T cells, the T cell clonality, the density of regulatory T cells
in the invasion margin, the alpha diversity of the microbiota and blood urea nitrogen.
This step resulted in a dataset with 336 variables and 26% of missing values.

Within tissue biomarkers, clusters of correlation were consistent with the hierarchical
structure of the cell populations (e.g., intra-tumor sub-populations of NK or myeloid-
derived suppressor cells) and location (e.g., T lymphocytes in the tumor, parenchyma
or stroma), see website. Within blood biomarkers, similar clusters were present with,
e.g., one cluster of lymphocytes subdivided into two clusters for T and B cells. No major
correlations were observed between circulating and tumor biomarkers, which clustered
mostly together. To mitigate multicollinearity for the subsequent biostatistical and
machine learning analyses, we used iterative variance inflation factor analysis⁴⁷. This
resulted in a decorrelated dataset with 214 features (Supplementary Table 1,
[ website ).

The primary outcome of this study was PrR to anti-PD(L)1 ICI⁴⁰,⁴¹. This was treated as a
binarized outcome of primary resistant or non-primary resistant. Based on this
definition, PrR to treatment was observed in 39.5% overall and in 29.9% (n=80/268)
and 55.1% (n=92/167) patients in the first- and subsequent-line settings, respectively
(Table 1). The secondary outcome of this study was PFS, measured from time of study
inclusion to progression of disease or death. Patients who had not progressed were
censored at the time of last contact with the study team. At a median duration of follow-
up of 16.9 months (95% CI 15.3-18.8), the median PFS was 10.1 months (95% CI 9.0-
14.3) for patients treated in the first-line setting. At a median duration of follow-up of
25.5 months (95% CI 21.2 – 29), median PFS was 3.77 months (95% CI 2.77- 5.43) for
patients treated in the subsequent-line setting (Supplementary Figure 2).

Biostatistical analysis revealed limited associations of single features with PrR
and PFS (Figure 2, Table S1)

Including patients treated in all treatment lines, 13 markers were significantly
associated with PrR, PFS and OS (PD-L1-controlled logistic/Cox regressions, FDR-
adjusted p-value < 0.05, **Figure 2A-B**). They included Regulatory T-cells (Treg) in the
stroma, C-reactive protein (CRP), Chlorine, Albumin, PD-L1 TPS (continuous), Red
blood cells (RBC), VCAM, ICAM, ECOG performance status, ASAT, Platelets, Sodium
and Alkaline phosphatase (**Figure 2B-D**). All the markers associated with PrR were
also associated with PFS except NLR, which was nevertheless associated with OS and
had an FDR-adjusted p-value of 0.0653. Thirty-three biomarkers showed significant
differences between patients with shorter versus longer PFS (FDR adjusted p-value
cut-point of <0.05, **Figure 2C-D**). Among these, 13 were also associated with OS and
12 biomarkers were specifically associated to OS (**Figure 2A, Supplementary Figure
3**).

However, no single feature emerged as being strongly associated with either PFS or
PrR and the predictive performance of individual features was modest at best, with the
highest AUC for association with PrR observed for the stroma Treg (odds ratio (OR) for
PrR 0.572, 95% CI 0.444 - 0.737, AUC 0.643), VCAM (odds ratio (OR) for PrR 2.15, 95% CI
1.36-3.4, AUC 0.64) and haemoglobin (odds ratio (OR) for PrR 0.623, 95% CI 0.493 - 0.787,
AUC 0.635). The highest C-index (association with PFS) was observed for the following
features: CRP (hazard ratio (HR) for PFS 1.52 (95% CI 1.31-1.76, C-index 0.61), chlorine
(HR 0.685, 95% CI 0.593-0.792, C-index 0.61), haemoglobin (HR 0.733, 95% CI 0.652-0.824,
C-index 0.61), RBC (HR 0.673, 95% CI 0.596-0.774, C-index 0.61), and PD1+ cells (HR
0.72, 95% CI 0.602-0.862, C-index 0.61) (**Figure 2B**, website).

PD-L1 TPS is a standard approved clinical biomarker for front-line ICI therapy in
NSCLC.¹ Using a cut point of TPS <1% compared with 1%, odds of PrR was
significantly higher in tumors with TPS <1% [n=200/362, OR 2.03 (95% CI 1.23-3.34),
p=5.29e-03] on univariable analysis, however this difference was non-significant when
adjusting for FDR (p-adjusted = 0.0933). When PD-L1 TPS was evaluated as a
continuous variable, this feature remained significantly associated with non-PrR
[n=362, OR 0.672 (95% CI 0.533-0.848), p=7.95e-04, p-adjusted 0.0319]. However,
PD-L1 TPS showed only moderate ability to discriminate between groups in the PrR
outcome (AUC 0.62) (**Figure 2A-B**). Out of the large proportion of patients with PD-L1
TPS < 1% only 52% (n=104/200) were PrR, indicating that TPS is a poor predictor of
resistance. Conversely, 71.2% (n=116/163) of the patients with TPS ≥ 1% were non-
PrR, highlighting the predictive value of TPS for response. Similarly, for the secondary
outcome of PFS, TPS <1% compared with 1% was associated with a significantly
shorter PFS on univariable analysis [n=200/365, HR 1.53 (95% CI 1.12-2.09),
p=7.01e-03], however this was non-significant on multivariable analysis (p-adjusted
0.0561). Higher PD-L1 TPS expression as a continuous variable remained significantly
associated with improved PFS [n=365, HR 0.661 (95% CI 0.566-0.772), p=1.52e-07, p-
adjusted 1.12e-05]. Again, the C-index for PD-L1 was poorly discriminative (0.60) (Figure 2A-B).

We then assessed features derived from WES data obtained from tumor tissue with matched normal specimens (blood) to detect somatic gene alterations. On PD-L1-controlled analysis in all patients treated with ICI regardless of line of treatment, no single gene mutation was significantly associated with PrR and only mutations in one gene, *CTNNB1*, were significantly associated with shorter PFS [n=17/309, HR 2.49 (95% CI 1.44-4.31), p-value = 1.08e-03, p-adjusted = 0.0158] (Figure 2C). Significance of mutations on cytokeratin-associated genes in Figure 2C was due to only one mutated patient with a very rapid progression (36 days). Mutations in *STK11* or *KEAP1* have previously been identified as being associated with ICI resistance in advanced NSCLC<sup>29</sup>. However, in this study, mutation in *STK11* was neither significantly associated with PrR [n=58/305, OR 1.17 (95% CI 0.634-2.17), AUC 0.54, p-adjusted = 0.896], nor with PFS [n=309, HR 1.38 (95% CI 0.972-1.95), C-index 0.542, p-adjusted = 0.26]. Similarly, mutation in *KEAP1* was not significantly associated with PrR [n=53/305, OR 0.663 (0.344-1.28), p-adjusted 0.58], nor PFS [n=309, HR 1.06 (95% CI 0.742-1.52), p-adjusted = 0.914] (website). High TMB is another approved biomarker for the selection of patients for ICI therapy in both advanced NSCLC and pan-tumor indications<sup>41</sup>, however in this study, TMB was non-discriminatory for the outcomes of PrR [n=232, OR 1.07 (95% CI 0.826-1.4), AUC 0.55, p = 0.592], nor PFS [n=234, HR 1.08 (95% CI 0.931-1.25), C-index 0.538, p=0.308] (dashboard).

## Development of a signature and supervised machine learning model for the prediction of PrR and PFS (Figure 1D, Figure 3)

To identify a stable, minimal biomarker signature predictive of PrR to anti-PD-(L)1 ICIs, we implemented a two-stage supervised learning framework comprising feature selection and predictive modeling (Figure 1D). Data from first- and later-line treatments were combined to maximize sample size (n = 435), while keeping the line of treatment as a predictor. Thirty-six feature selection methods were benchmarked to identify a reduced subset predictive of PrR. These were combined with 10 machine learning models, resulting in 360 feature selection-machine learning (FS-ML) combinations<sup>46</sup>. Strong emphasis was set on the stability of the selected subset (i.e., the similarity between selected subsets obtained under small data perturbations),
assessed using a metric defined on 100 subsamples of the data⁴⁶,⁵². The standard LASSO method (and adaptations), as well as most machine learning-based feature selection algorithms, exhibited weak stability results (< 0.6, website).

The feature selection method that demonstrated the best stability (0.69) and top AUC for prediction of PrR was Benjamini-Hochberg FDR-adjusted statistical filtering of p-values from PFS Cox regression at the 0.01 level. It resulted in a refined set of 18 features, hereafter referred to as the “18-feature signature”. For predictive modeling, a gradient boosting model achieved the best .632 optimism-corrected AUC (0.73, website), and second-best .632 optimism-corrected PPV at 60.5%, after an extreme gradient boosting at 63.2%, but with much lower bias (difference between apparent and corrected PPV) of 0.176 (versus 0.35). For further details of the machine learning system, see Methods.

## Supervised machine learning for prediction of PrR (Figure 3)

Two supervised machine learning models were developed to predict PrR to anti-PD-(L)1 therapy. The first model was trained on the broad set of 214 decorrelated features with no further feature selection, while the second was trained on the 18 features signature (Figure 3A). Both models generated continuous risk scores ranging from 0 to 1, with higher scores reflecting increased predicted probability of PrR. The primary performance metrics were positive predictive value (PPV) and AUC. Secondary metrics included balanced accuracy (ACC), sensitivity (SE), specificity (SP), and negative predictive value (NPV). All metrics were assessed using optimism-correction over 100 subsampled datasets. The 214-feature model achieved a PPV of 63%, AUC of 0.73 and ACC of 68%, outperforming any individual feature previously assessed in univariable analyses (Figure 3A). PPV is to be appreciated compared to the 39.5% proportion of PrR patients (red bar in Figure 3A) and PPVs of 49.6% for PD-L1 (threshold 1%) and 43.1% TMB (threshold 10 mutations/megabases), respectively. The 18-feature signature demonstrated comparable performance, with a more favorable sensitivity-specificity trade-off, suggesting improved generalizability and interpretability. However, PPV remained the limiting metric (0.63 and 0.60 for the 214- and 18-feature models, respectively), underscoring the intrinsic difficulty of accurately predicting PrR using static baseline data.
The 18-feature signature was selected for further validation due to its parsimony and robust performance. When stratified by treatment line, the model achieved PPVs of 51% and 64% for first-line and subsequent-line ICI-treated subgroups, respectively (to be compared to 29.8% and 55.1% of PrR patients in the two populations, respectively). In first-line patients, secondary model performance metrics were: AUC 69%, NPV 79%, SE 47%, and SP 81% (Figure 3A). In subsequent-line patients, the signature achieved AUC 69%, NPV 63%, SE 77%, and SP 45% (Figure 3A).

Notably, although the model was trained solely to predict binary PrR outcomes, it also demonstrated strong prognostic value for PFS, supporting its clinical relevance (Figure 3B). Patients were dichotomized into high-risk and low-risk groups based on the model-derived score. In the training cohort (n = 289), high-risk patients had a median PFS of 2.37 months (95% CI: 1.47–2.8) compared to 18.6 months (95% CI: 13.1–26.3) in the low-risk group (HR = 0.0904, 95% CI: 0.0642–0.127; log-rank p < 0.0001). Similar stratification was observed in the test set (n = 145), where high-risk patients had a median PFS of 3.9 months (95% CI: 1.63–5.27) versus 14.6 months (95% CI: 9.77–NA) for low-risk patients (HR = 0.307, 95% CI: 0.204–0.462; p < 0.0001) (Figure 3C). Within the test set, subgroup analyses further supported model generalizability. Among first-line treated patients (n = 88), high-risk individuals had a median PFS of 4.53 months (95% CI: 2.43–NA) versus 14.6 months (95% CI: 9.77–NA) in the low-risk group (HR = 0.392, 95% CI: 0.202–0.761; p = 0.0042). In subsequent-line patients, high-risk individuals had a median PFS of 2.77 months (95% CI: 1.4–5.27), compared to 14.1 months (95% CI: 6.33–NA) in the low-risk group (HR = 0.315, 95% CI: 0.157–0.63; p = 0.00062) (Figure 3B). These results highlight the robustness and prognostic potential of the 18-feature signature, reinforcing its utility for identifying patients at high risk of early disease progression following ICI treatment.

The predictive capacity of the final 18-feature multivariable machine learning signature was benchmarked against leading univariable biomarkers, including PD-L1 TPS and TMB. The machine learning signature consistently outperformed current standard-of-care biomarkers. Specifically, in the first-line setting, PD-L1 TPS using a cut-point of 1% yielded a PPV of 36% (33% with cut-point at 50%) and PD-L1 continuous a .632 optimism-corrected AUC of 0.56 (Figure 3C), comparable to previously reported performance metrics (AUC ≈ 0.60¹⁷; PPV ≈ 40%¹⁸). Similarly, TMB demonstrated lower discriminative ability than the machine learning signature, with TMB (cut-off 10
megabases) yielding a PPV of 35% and an AUC of 58% in the first-line setting. In comparison, PPV of the machine learning signature was 51% and AUC 69% (Figure 3A-C). These findings underscore the limited predictive utility of individual biomarkers in isolation and highlight the added value of integrative, multivariable modeling in capturing the complex biological underpinnings of ICI resistance.

To assess the added value of the integrated machine learning signature, we compared its predictive performance to that of individual biomarkers grouped by type (tumor biomarkers, circulating immune cells, vasculo-monitoring, routine laboratory blood tests). For each group, performances were obtained applying the FS + ML pipeline described above. When benchmarked against grouped biomarker types the machine learning signature demonstrated significantly superior performance across multiple evaluation metrics, particularly PPV, AUC and ACC (Figure 3D), underscoring the limitations of unimodal predictors and the benefit of data integration. Among the individual biomarker data types, routine blood test parameters emerged as the most informative contributors to model performance, highlighting their underutilized potential in predictive immuno-oncology. These findings reinforce the hypothesis that orthogonal, multi-modal data sources capture complementary aspects of tumour-host biology, and that their integration through machine learning yields more accurate and clinically actionable biomarkers than any single modality alone.

## Multivariable signature for prediction of PrR (Figure 4)

The final machine learning signature comprised 18 features (excluding the line of treatment) spanning multiple biomarker types, including genomic, clinical, immunologic, and laboratory-based parameters (Figure 4A-B). Feature contribution analysis using SHapley Additive exPlanations (SHAP) revealed complex nonlinear relationships between individual features and model output, reflecting the multidimensional biological underpinnings of PrR to ICI. Among the top-ranked predictors, the most discriminative feature after PD-L1 was stromal regulatory T cell (Treg) infiltration, with low Treg abundance strongly associated with resistance. Despite efforts to reduce multicollinearity and apply stringent variable selection, residual correlations persisted within the signature, particularly between biomarkers that are classically associated (e.g., chloride and sodium, r = 0.663 or albumin and CRP, r = -0.515, Pearson correlation, website). This indicates that, despite their
correlation, these variables still contributed complementary predictive information when combined (Supplementary Figure 4).

Strikingly, variables from routine blood tests accounted for approximately half of the final signature (Figure 4A-B). Chloride concentration emerged as the most discriminating laboratory variable and was the top split in the model-derived decision tree (Figure 4C), followed by albumin, RBC, monocytes-to-lymphocytes ratio (MLR), alkaline phosphatase and CRP—features also independently validated in prior studies.³⁵ PD-L1 expression was an important contributor to the model but exhibited a non-monotonic, nonlinear effect (Figure 4A), possibly explaining why it did not appear in the upper nodes of the decision trees (Figure 4C).

## A machine learning-based digital tool to predict individual risk of resistance (Figure 5)

To interrogate the individual contributions of each biomarker within the multivariable machine learning signature, we further employed SHAP analysis at the individual scale to provide a digital tool that not only predicts individual risk of resistance but also explains which variables contribute most to favoring or preventing resistance. The model was the one fitted on the training dataset and predictions were made for patients in the test set (Figure 5). This approach allowed us to disentangle the distinct and often nonlinear influence of individual biomarkers, even within a highly interdependent multivariable model. For instance, in some patients, low stromal Treg infiltration was a primary risk driver, preceding low PD-L1 TPS (Figure 5A-B). Notably, while some signature factors such as platelets or chlorine differed strongly between these two patients, this did not significantly alter prediction risks. In other patients, systemic inflammation parameters emerged as primary risk drivers, such as high CRP (Figure 5C). In this PrR patient, serum CRP was markedly elevated (CRP = 235 mg/L), and his resistance risk was predicted to be 77.7%, despite high PD-L1 TPS (100%) and high stromal Treg infiltration (164 cells/mm²). In another PrR patient example (Figure 5D) – correctly predicted by the model with a resistance risk of 74.5% – elevated MLR (MLR = 1.5) was among the main resistance drivers, despite purportedly protective factors such as low stroma Treg and high CD3+CD4+ infiltration. Some long-term responsive patients could be explained by high PD-L1 TPS (Figure 5E), but this was not always the main contributive factor, compared with, e.g., nutritional status quantified by albumin (Figure 5F). On the other hand, some patients had low PD-L1 TPS (PD-L1
expression = 5%) but still exhibited ≥ 12 months response (Figure 5G). The machine
learning model explained this by a relatively high level of alkaline phosphatase,
chlorine and ASAT, in the addition to a favorable MLR and non-zero TPS. Another
patient exhibiting long term stable disease, despite being in the subsequent-line group,
was found to be consistently predicted non-PrR, with one of the main drivers identified
as protective of progression being the level of sodium (Figure 5H).

These individualized SHAP profiles provided critical insight into the
heterogeneity of resistance mechanisms and illustrated how specific biomarker
contributions varied not only in magnitude but also in direction across patients. Such
interpretability is essential for both biological understanding and future clinical
translation of machine learning -derived predictive models.
DISCUSSION

Despite extensive research into host-tumour interactions in the immune regulation of advanced NSCLC, confirmed biomarkers predictive of primary ICI resistance remain elusive. This study aimed to investigate a broad spectrum of potential biomarkers by performing a comprehensive assessment of 443 clinical factors, routine laboratory blood tests, tumour microenvironmental characteristics and blood-derived immune and vascular markers in a cohort of 439 advanced NSCLC patients prior to the receipt of anti-PD(L)1 ICI. We released here an interactive dashboard that enables clinicians and researchers to explore this extensive data.

The integration of multi-parametric biomarkers to guide oncologic decision-making remains an emerging field. In this study, we demonstrated that a machine learning–based, multimodal signature combining data cited above can improve prediction of PrR to anti-PD-(L)1 ICI in patients with advanced NSCLC compared with standard unimodal biomarkers. Some but few prior studies have similarly utilized machine learning approaches to integrate clinical, radiomic, pathologic and genomic data for ICI outcome prediction, combining biomarkers within or across modalities <sup>19,35–38,77</sup>. In comparison, our work brings novel contributions in several aspects.

First, we addressed a biologically complex and clinically urgent question by focusing on PrR, following the SITC guidelines <sup>38,39</sup>. Previous studies typically use best overall response (BOR) or OS as a primary endpoints<sup>36,37</sup> — metrics that depend on follow-up duration and are generally less challenging to predict. Resistance to therapy is inherently more difficult to predict (e.g., PD-L1 TPS ≥ 50% already serves as a reasonable predictor of response), and studies addressing PrR have reported only AUC values with modest performance, and not clinically actionable metrics such as PPV for resistance<sup>35,38</sup>.

Second, our study included the largest cohort of patients with multi-modal data reported to date, with a substantial proportion treated with combination therapy – the current standard of care in first-line advanced NSCLC. Previous studies with significantly larger sample sizes have focused on single modalities (e.g., pathology<sup>37</sup> or clinical<sup>35,78</sup>) or included datasets with substantially fewer patients<sup>36,38</sup>.

Third, we devoted rigorous attention to feature selection (FS), aiming to identify a sparse and stable signature predictive of PrR. This challenge is frequently overlooked
in the biomedical literature, where multi-modal models often integrate all available
variables or rely on a single feature selection method applied only to the training set or
entire dataset³⁷,³⁸,⁷⁸. In contrast, we systematically benchmarked 36 FS methods for
both stability and performance, using an optimism-correction framework that controls
overfitting bias⁵¹,⁵².

The 18-feature machine learning model integrated a diverse array of biomarker types, including histologic, clinical, immunologic, and biochemical parameters. Notably, half of the signature derived from this model comprises blood-based parameters readily available in routine clinical practice. Serum chloride, albumin, MLR, alkaline phosphatase, and CRP emerged among the top discriminators, echoing recent reports validating the prognostic value of routine markers in immuno-oncology³⁵,⁷⁸ and extend these findings by situating such features within a robust, multivariable predictive framework. Overall, the prominence of blood-based biomarker features that were selected within the machine learning model may reflect their ability to capture systemic immune and metabolic changes, and offer a potential advantage over tissue-based markers, which are often limited by sampling bias and the spatial tumour heterogeneity.

Our data reaffirm the limitations of existing, tissue-centric biomarkers: PD-L1 expression and TMB demonstrated only modest discrimination (AUCs ~0.60, PPVs ~35% in the first line setting), and canonical resistance mutations in *STK11* or *KEAP1* failed to retain significance after correction for multiple testing. In contrast, the ML signature achieved an optimism-corrected AUC of 0.69 and positive predictive values of 51% and 64% in first- and later-line settings, respectively — substantially improving upon baseline PrR rates of 29.9% and 55.1%. Encouragingly, the multimodal model robustly stratified early non-responders, as evidenced by sharp separation of Kaplan–Meier progression-free survival curves early after treatment initiation, in the training and tests sets. This suggests that our signature holds potential utility in the early identification of patients most likely to be primary resistant to anti-PD(L)1 ICI, and who may benefit from treatment intensification approaches, closer clinical monitoring for disease progression or exploration of alternative combinations (e.g., antibody-drug conjugates), which should be prospectively validated in future clinical trials.
Beyond prediction, the interpretability of the model enabled by tools such as SHAP
values offered a unique window into the underlying biology of ICI resistance. By
quantifying the patient-specific contributions of immune, metabolic, and clinical
features, the model can generate mechanistically plausible hypotheses about causal
pathways and therapeutic vulnerabilities. However, the absence of non-ICI-treated
patients prevented us to distinguish between prognostic and predictive value of the
biomarkers, either alone or combined. Nevertheless, the signature may also serve as a
foundation for future work in causal inference and systems-level understanding of
immunotherapy resistance.

While the present analysis was limited to biomarker data obtained at pre-treatment
timepoints, the FS + ML pipeline described here offers a scalable framework for future
analyses of longitudinally collected biomarker datasets. Applying this approach to
samples obtained during the course of first-line ICI therapy could provide critical
insights into the dynamic mechanisms that underlie not only primary resistance but also
the development of acquired resistance. Emerging evidence suggests that therapeutic
pressure during first-line treatment can profoundly reshape the tumor-immune
ecosystem<sup>79</sup>. This includes the expansion of immunosuppressive cell subsets<sup>79</sup>, clonal
evolution of antigen-deficient tumor populations<sup>80</sup>, and the induction of compensatory
immune checkpoints that contribute to T cell exhaustion<sup>30</sup>. Concomitant systemic
alterations, such as lymphodepletion and inflammatory cytokine shifts, further dampen
immune-mediated tumor control<sup>79</sup>. These evolving biological features underscore the
importance of real-time biomarker monitoring and support the development of adaptive
treatment strategies designed to intercept resistance mechanisms as they emerge.

Limitations of this study include its restriction to baseline data, which limits insight into
the temporal evolution of predictive features. Furthermore, while the optimism-
corrected performance metrics and internal resample validation reduce the risk of
overfitting, the lack of external validation in an independent cohort necessitates
cautious interpretation. Additionally, generalizability to external populations may be
constrained by the inclusion of exploratory biomarkers, such as tumour
microenvironmental or circulating vascular features, that are not routinely captured in
standard care. An additional limitation of our study remains that the integration of
heterogeneous data sources requires non-trivial infrastructure and technical expertise
to harmonise disparate data formats, ensure data quality and maintain consistent
annotation across platforms. While these technical demands certainly may act as a
barrier to implementation in real-world clinical settings, they also reflect the scale of
nation-wide collaboration that was required to complete this ambitious project.

Despite these challenges, our findings provide proof-of-principle that multimodal
signatures derived from myriad clinical input and translational assays can meaningfully
augment ICI resistance prognostication. To facilitate broader adoption, we provide the
underlying workflows and software pipelines for use within the academic oncology
community. These tools are adaptable beyond NSCLC and may enable scalable
deployment of similar integrative models in other cancer types and disease settings.
Future directions include prospective validation of the 18-feature ML signature in
external cohorts, incorporation of longitudinal data to capture resistance dynamics over
time, and integration of emerging data types such as microbiome profiles. Ultimately,
embedding such predictive frameworks into future clinical trials could enable
biomarker-guided immunotherapy selection, personalized risk-adapted strategies, and
more rational trial design for overcoming resistance in NSCLC.
# FIGURE LEGENDS

## Figure 1: Workflow of the analysis

A. Multi-modal biomarkers

B. Quantification of the number of biomarkers by modality

C. Primary endpoint: primary resistance

D. Schematic of the machine learning feature selection and predictive modeling. NA = not available (missing values).

## Figure 2: Univariable biostatistics

A. Venn diagram of the significant markers for primary resistance (PrR), PFS and OS

B. Univariable statistics of markers significantly associated with PrR, PFS and OS.

C. Forest plot of markers significantly associated with PFS (PD-L1-controlled Cox regression, FDR-adjusted p-value < 0.05). Uncorrelated dataset.

D. Box and Kaplan-Meier (KM) plots of the most significant features for PFS (log-rank test), for each type of biomarkers. Red dashed line = KM optimal cutoff for the Kaplan-Meier.

Unless expressed differently, all statistics are from PD-L1-controlled logistic/Cox regressions, FDR-adjusted p-value < 0.05 on the uncorrelated dataset.

See the dashboard for direct exploration of all individual biomarkers at https://compo.inria.fr/pioneer-dashboard/.

## Figure 3: Machine learning predictive performances

A. Predictive metrics of PrR, trained on all the uncorrelated biomarkers (left) or the 18-biomarker signature (right). Bars = standard deviations across subsamples.

B. Progression-free survival stratification performances of the signature + machine learning model on a training (left) and test dataset. The cut point was defined on the training dataset and propagated to the test dataset.
C. Optimism-corrected metrics for the signature + ML model versus PD-L1 expression (1% threshold for all metrics except AUC where the continuous marker was used) and tumor mutational burden (10 megabases threshold for all metrics except AUC where the continuous marker was used).

D. Comparative metrics of PrR between PD-L1, each BM type and the signature. Bars = standard deviations across subsamples.

1L = first-line patients (combination therapy); 2+L = subsequent-line patients (monotherapy); ACC = balanced accuracy; AUC = area under the receiver-operator curve; NPV = negative predictive value; PPV = positive predictive value; SE = sensitivity; SP = specificity. Positive event = primary resistance (PrR).

**Figure 4: Machine learning signature**

A. Signature resulting from feature selection on the decorrelated biomarkers dataset, sorted by random forest SHAP global importance. One point is the SHAP value of one patient for one variable of the signature. The variable of the line of treatment was not considered as a biomarker and is not reported here.

B. Hierarchical clustering performed on the biomarkers dataset, grouped by line setting (first line: combination therapy, second+ line: monotherapy) and primary resistance.

C. Decision trees using signature.

MLR: monocytes-to-lymphocytes ratio. VCAM: Vascular Cell Adhesion Molecule. ICAM: Inter-Cellular Adhesion Molecule. ASAT: Aspartate Aminotransferase. ECOG: Eastern Cooperative Oncology Group performance status.

**Figure 5: Digital tool to predict individual risk of resistance**

Left: Individual SHAP waterfall plots representing the contributions of each biomarker from to the final prediction. Departing from a basal risk in the population (reflecting the proportion of PrR patients) indicated as E[f(x)], each feature contributes either positively (red) or negatively to the risk of PrR, leading to an individual risk denoted f(x). The model was fitted on the training dataset and predictions shown here were performed on patients from the test set. Middle radar plots show scaled values of the features of the patient, expressed as percentiles in the total population. Right plots
show time evolution of the patient SLD, with indications of RECIST evaluations and
progression / death when relevant.

A-B : PrR patients

C-D : non-PrR patients
ACKNOWLEDGEMENTS

This work benefited from a government grant handled by the French National Research Agency (ANR) as part of the France 2030 investment plan, under the reference ANR-17-RHUS-0007.

This work was supported by a partnership of Aix-Marseille Université (AMU), Assistance Publique Hôpitaux de Marseille (APHM), Centre National de La Recherche Scientifique (CNRS), Institut National de la Santé et de la Recherche Médicale (INSERM), Centre Léon Bérard (CLB), Institut Paoli Calmettes (IPC), Gustave Roussy (GR), AstraZeneca (AZ), Veracyte (VERA), Innate Pharma (IPH) & ImCheck Therapeutics (ICT), and initiated by Marseille Immunopole.

The authors gratefully acknowledge the support of the APHM, which sponsored the PIONeeR clinical studies. Its role was to control the appropriateness of ethical and legal considerations for all centers and to perform the monitoring of the consents signed and the clinical data recorded and coded as part of the study.

The authors are grateful to all the patients and their families, as well as all the investigators, for their participation in the study.

This work benefited from support from ITMO Cancer AVIESAN and French Institut National du Cancer (grant #19CM148-00) and from the French National Research Agency (ANR), under the France 2030 program, reference ANR-22-PESN-0017.
REFERENCES

1. Hendriks, L. E. et al. Non-oncogene-addicted metastatic non-small-cell lung cancer: ESMO Clinical Practice Guideline for diagnosis, treatment and follow-up. *Ann Oncol* **34**, 358–376 (2023).

2. Wang, Y. et al. Immunotherapy for advanced-stage squamous cell lung cancer: the state of the art and outstanding questions. *Nat Rev Clin Oncol* **22**, 200–214 (2025).

3. Gandhi, L. et al. Pembrolizumab plus Chemotherapy in Metastatic Non-Small-Cell Lung Cancer. *N. Engl. J. Med.* **378**, 2078–2092 (2018).

4. Herbst, R. S. et al. Pembrolizumab versus docetaxel for previously treated, PD-L1-positive, advanced non-small-cell lung cancer (KEYNOTE-010): a randomised controlled trial. *The Lancet* **387**, 1540–1550 (2016).

5. Mazieres, J. et al. Atezolizumab Versus Docetaxel in Pretreated Patients With NSCLC: Final Results From the Randomized Phase 2 POPLAR and Phase 3 OAK Clinical Trials. *J Thorac Oncol* **16**, 140–150 (2021).

6. Reck, M. et al. Pembrolizumab versus Chemotherapy for PD-L1-Positive Non-Small-Cell Lung Cancer. *N. Engl. J. Med.* **375**, 1823–1833 (2016).

7. Rittmeyer, A. et al. Atezolizumab versus docetaxel in patients with previously treated non-small-cell lung cancer (OAK): a phase 3, open-label, multicentre randomised controlled trial. *The Lancet* **389**, 255–265 (2017).

8. Socinski, M. A. et al. Atezolizumab for First-Line Treatment of Metastatic Nonsquamous NSCLC. *N. Engl. J. Med.* **378**, 2288–2301 (2018).

9. Borghaei, H. et al. Five-Year Outcomes From the Randomized, Phase III Trials CheckMate 017 and 057: Nivolumab Versus Docetaxel in Previously Treated Non-Small-Cell Lung Cancer. *J Clin Oncol* **39**, 723–733 (2021).
10. Mok, T. S. K. et al. Associations of tissue tumor mutational burden and mutational status with clinical outcomes in KEYNOTE-042: pembrolizumab versus chemotherapy for advanced PD-L1-positive NSCLC. *Ann Oncol* **34**, 377–388 (2023).

11. Mok, T. S. K. et al. Pembrolizumab versus chemotherapy for previously untreated, PD-L1-expressing, locally advanced or metastatic non-small-cell lung cancer (KEYNOTE-042): a randomised, open-label, controlled, phase 3 trial. *Lancet* **393**, 1819–1830 (2019).

12. Reck, M. et al. Nivolumab plus ipilimumab versus chemotherapy as first-line treatment in advanced non-small-cell lung cancer with high tumour mutational burden: patient-reported outcomes results from the randomised, open-label, phase III CheckMate 227 trial. *Eur J Cancer* **116**, 137–147 (2019).

13. Ricciuti, B. et al. Association of High Tumor Mutation Burden in Non-Small Cell Lung Cancers With Increased Immune Infiltration and Improved Clinical Outcomes of PD-L1 Blockade Across PD-L1 Expression Levels. *JAMA Oncol* **8**, 1160–1168 (2022).

14. Cristescu, R. et al. Pan-tumor genomic biomarkers for PD-1 checkpoint blockade–based immunotherapy. *Science* **362**, eaar3593 (2018).

15. Yarchoan, M., Hopkins, A. & Jaffee, E. M. Tumor Mutational Burden and Response Rate to PD-1 Inhibition. *N Engl J Med* **377**, 2500–2501 (2017).

16. Hellmann, M. D. et al. Genomic Features of Response to Combination Immunotherapy in Patients with Advanced Non-Small-Cell Lung Cancer. *Cancer Cell* **33**, 843-852 e4 (2018).

17. Rizvi, H. et al. Molecular Determinants of Response to Anti-Programmed Cell Death (PD)-1 and Anti-Programmed Death-Ligand 1 (PD-L1) Blockade in Patients
With Non-Small-Cell Lung Cancer Profiled With Targeted Next-Generation Sequencing. *J. Clin. Oncol.* https://doi.org/10.1200/JCO.2017.75.3384 (2018) doi:10.1200/JCO.2017.75.3384.

18. Camidge, D. R., Doebele, R. C. & Kerr, K. M. Comparing and contrasting predictive biomarkers for immunotherapy and targeted therapy of NSCLC. *Nat. Rev. Clin. Oncol.* **16**, 341–355 (2019).

19. Chowell, D. et al. Patient HLA class I genotype influences cancer response to checkpoint blockade immunotherapy. *Science* **359**, 582–587 (2018).

20. Postow, M. A. et al. Peripheral T cell receptor diversity is associated with clinical outcomes following ipilimumab treatment in metastatic melanoma. *J Immunother Cancer* **3**, 23 (2015).

21. Reuben, A. et al. Comprehensive T cell repertoire characterization of non-small cell lung cancer. *Nat Commun* **11**, 603 (2020).

22. Sun, R. et al. A radiomics approach to assess tumour-infiltrating CD8 cells and response to anti-PD-1 or anti-PD-L1 immunotherapy: an imaging biomarker, retrospective multicohort study. *Lancet Oncol.* **19**, 1180–1191 (2018).

23. Lee, J. S. & Ruppin, E. Multiomics Prediction of Response Rates to Therapies to Inhibit Programmed Cell Death 1 and Programmed Cell Death 1 Ligand 1. *JAMA Oncol* **5**, 1614–1618 (2019).

24. Higgs, B. W. et al. Interferon Gamma Messenger RNA Signature in Tumor Biopsies Predicts Outcomes in Patients with Non-Small Cell Lung Carcinoma or Urothelial Cancer Treated with Durvalumab. *Clin Cancer Res* **24**, 3857–3866 (2018).

25. Ayers, M. et al. IFN-gamma-related mRNA profile predicts clinical response to PD-1 blockade. *J Clin Invest* **127**, 2930–2940 (2017).
26. Vanhersecke, L. et al. Mature tertiary lymphoid structures predict immune checkpoint inhibitor efficacy in solid tumors independently of PD-L1 expression. *Nat Cancer* **2**, 794–802 (2021).

27. Sade-Feldman, M. et al. Resistance to checkpoint blockade therapy through inactivation of antigen presentation. *Nat Commun* **8**, 1136 (2017).

28. Shin, D. S. et al. Primary Resistance to PD-1 Blockade Mediated by JAK1/2 Mutations. *Cancer Discov* **7**, 188–201 (2017).

29. Skoulidis, F. et al. Co-occurring genomic alterations define major subsets of KRAS-mutant lung adenocarcinoma with distinct biology, immune profiles, and therapeutic vulnerabilities. *Cancer Discov* **5**, 860–77 (2015).

30. Koyama, S. et al. STK11/LKB1 Deficiency Promotes Neutrophil Recruitment and Proinflammatory Cytokine Production to Suppress T-cell Activity in the Lung Tumor Microenvironment. *Cancer Res* **76**, 999–1008 (2016).

31. Kitajima, S. et al. Suppression of STING Associated with LKB1 Loss in KRAS-Driven Lung Cancer. *Cancer Discov* **9**, 34–45 (2019).

32. Han, G. et al. 9p21 loss confers a cold tumor immune microenvironment and primary resistance to immune checkpoint therapy. *Nat Commun* **12**, 5606 (2021).

33. Bessede, A. et al. TROP2 Is Associated with Primary Resistance to Immune Checkpoint Inhibition in Patients with Advanced Non-Small Cell Lung Cancer. *Clin Cancer Res* **30**, 779–785 (2024).

34. Kalbasi, A. & Ribas, A. Tumour-intrinsic resistance to immune checkpoint blockade. *Nat Rev Immunol* **20**, 25–39 (2020).

35. Yoo, S.-K. et al. Prediction of checkpoint inhibitor immunotherapy efficacy for cancer using routine blood tests and clinical data. *Nat. Med.* **31**, 869–880 (2025).
36. Vanguri, R. S. et al. Multimodal integration of radiology, pathology and genomics for prediction of response to PD-(L)1 blockade in patients with non-small cell lung cancer. *Nat Cancer* **3**, 1151–1164 (2022).

37. Rakaee, M. et al. Deep Learning Model for Predicting Immunotherapy Response in Advanced Non-Small Cell Lung Cancer. *JAMA Oncol* **11**, 109–118 (2025).

38. Captier, N. et al. Integration of clinical, pathological, radiological, and transcriptomic data improves prediction for first-line immunotherapy outcome in metastatic non-small cell lung cancer. *Nat. Commun.* **16**, 614 (2025).

39. Eisenhauer, E. A. et al. New response evaluation criteria in solid tumours: Revised RECIST guideline (version 1.1). *Eur. J. Cancer* **45**, 228–247 (2009).

40. Rizvi, N. et al. Society for Immunotherapy of Cancer (SITC) consensus definitions for resistance to combinations of immune checkpoint inhibitors with chemotherapy. *J Immunother Cancer* **11**, (2023).

41. Atkins, M. B. et al. Society for Immunotherapy of Cancer (SITC) consensus definitions for resistance to combinations of immune checkpoint inhibitors with targeted therapies. *J Immunother Cancer* **11**, (2023).

42. Vallier, L. et al. Increasing the sensitivity of the human microvesicle tissue factor activity assay. *Thromb Res* **182**, 64–74 (2019).

43. Chalmers, Z. R. et al. Analysis of 100,000 human cancer genomes reveals the landscape of tumor mutational burden. *Genome Med* **9**, 34 (2017).

44. Nguyen Phuong, L., Vaglio, A. & Benzekry, S. compo.EDA: Data Analysis for Oncology. https://gitlab.inria.fr/compo/compo.eda. (2025).

45. Vaglio, A. & Benzekry, S. compo.tidyML. https://gitlab.inria.fr/compo/compo.tidyml. (2025).
46. Bakhmach, A. et al. ROOFS: RObust biOmarker Feature Selection. Preprint at
https://doi.org/10.48550/arXiv.2601.05151 (2026).

47. Gareth, G., Witten, D., Hastie, T. & Tibshirani, R. *An Introduction to Statistical Learning: With Applications in R*. (Springer US, 2021).

48. Yekutieli, D. & Benjamini, Y. Resampling-based false discovery rate controlling multiple test procedures for correlated test statistics. *J. Stat. Plan. Inference* **82**, 171–196 (1999).

49. Pedregosa, F. et al. Scikit-learn: Machine Learning in Python. *J Mach Learn Res* **12**, 2825–2830 (2011).

50. Efron, B. & Gong, G. A Leisurely Look at the Bootstrap, the Jackknife, and Cross-Validation. *Am. Stat.* **37**, 36–48 (1983).

51. Collins, G. S. et al. Evaluation of clinical prediction models (part 1): from development to external validation. *BMJ* **384**, e074819 (2024).

52. Nogueira, S. On the Stability of Feature Selection Algorithms. *J. Mach. Learn. Res.* **18**, 1–54 (2018).

53. Tibshirani, R. Regression Shrinkage and Selection via the Lasso. *J. R. Stat. Soc. Ser. B Methodol.* **58**, 267–288 (1996).

54. Meinshausen, N. & Bühlmann, P. Stability selection. *J. R. Stat. Soc. Ser. B Stat. Methodol.* **72**, 417–473 (2010).

55. Bach, F. R. Bolasso: model consistent Lasso estimation through the bootstrap. in *Proceedings of the 25th international conference on Machine learning* 33–40 (Association for Computing Machinery, New York, NY, USA, 2008).
doi:10.1145/1390156.1390161.

56. Zou, H. The Adaptive Lasso and Its Oracle Properties. *J. Am. Stat. Assoc.* **101**, 1418–1429 (2006).
57. Meinshausen, N. Relaxed Lasso. *Comput. Stat. Data Anal.* **52**, 374–393 (2007).

58. Simon, N., Friedman, J., Hastie, T. & Tibshirani, R. A Sparse-Group Lasso. *J. Comput. Graph. Stat.* **22**, 231–245 (2013).

59. Hédou, J. et al. Discovery of sparse, reliable omic biomarkers with Stabl. *Nat. Biotechnol.* **42**, 1581–1593 (2024).

60. Jenul, A. et al. RENT -- Repeated Elastic Net Technique for Feature Selection. *IEEE Access Vol 9 Pp 152333-152346 2021* **9**, 152333–152346 (2020).

61. Li, J. et al. Feature Selection: A Data Perspective. *ACM Comput. Surv.* **50**, 1–45 (2017).

62. Nie, F., Xiang, S., Jia, Y., Zhang, C. & Yan, S. Trace Ratio Criterion for Feature Selection. in *inproceedings* (eds Fox, D. & Gomes, C. P.) 671–676 ({AAAI} Press, 2008).

63. Lin, D. & Tang, X. Conditional Infomax Learning: An Integrated Framework for Feature Extraction and Fusion. in (eds Leonardis, A., Bischof, H. & Pinz, A.) 68–82 (Springer Berlin Heidelberg, 2006).

64. Fleuret, F. Fast Binary Feature Selection with Conditional Mutual Information. *J Mach Learn Res* **5**, 1531–1555 (2004).

65. Meyer, P. E. & Bontempi, G. On the Use of Variable Complementarity for Feature Selection in Cancer Classification. in (eds Rothlauf, F. et al.) 91–102 (Springer Berlin Heidelberg, 2006).

66. Jakulin, A. Machine learning based on attribute interactions. (Univerza v Ljubljani, 2005).

67. Yang, H. Feature selection based on joint mutual information. in vol. 23 (1999).

68. Bennasar, M., Hicks, Y. & Setchi, R. Feature selection using Joint Mutual Information Maximisation. *Expert Syst. Appl.* **42**, 8520–8532 (2015).
69. N. Spolaor, E. A. Cherman, M. C. Monard & H. D. Lee. ReliefF for Multi-label Feature Selection. in 6–11 (2013). doi:10.1109/BRACIS.2013.10.

70. Ding, C. & Peng, H. Minimum redundancy feature selection from microarray gene expression data. *J Bioinform Comput Biol* **3**, 185–205 (2005).

71. Xia, S. & Yang, Y. A model-free feature selection technique of feature screening and random forest based recursive feature elimination. (2023).

72. Calzolari, M. Shapicant. (2020).

73. Kuhn, M. & Wickham, H. *Tidymodels: A Collection of Packages for Modeling and Machine Learning Using Tidyverse Principles*. (2020).

74. Mayer, M. & Watson, D. kernelshap: Kernel SHAP. (2022).

75. Mayer, M. shapviz: SHAP Visualizations. (2022).

76. Lundberg, S. M. & Lee, S.-I. A unified approach to interpreting model predictions. in *Proceedings of the 31st International Conference on Neural Information Processing Systems* 4768–4777 (Curran Associates Inc., Red Hook, NY, USA, 2017).

77. Avallone, A. et al. Critical role of bevacizumab scheduling in combination with pre-surgical chemo-radiotherapy in MRI-defined high-risk locally advanced rectal cancer: results of the branch trial. *Oncotarget* **6**, 30394–30407 (2015).

78. Chang, T.-G. et al. LORIS robustly predicts patient outcomes with immune checkpoint blockade therapy using common clinical, pathologic and genomic features. *Nat. Cancer* **5**, 1158–1175 (2024).

79. Sharma, P., Hu-Lieskovan, S., Wargo, J. A. & Ribas, A. Primary, Adaptive, and Acquired Resistance to Cancer Immunotherapy. *Cell* **168**, 707–723 (2017).
80. Anagnostou, V. *et al.* Evolution of Neoantigen Landscape during Immune Checkpoint Blockade in Non-Small Cell Lung Cancer. *Cancer Discov* **7**, 264–276 (2017).
<table><thead><tr><th rowspan="2">Variable</th><th rowspan="2">Overall, N = 435<sup>1</sup></th><th colspan="2">Primary resistance</th><th rowspan="2">p-value<sup>3</sup></th></tr><tr><th>0, N = 263<sup>2</sup></th><th>1, N = 172<sup>2</sup></th></tr></thead><tbody><tr><td><strong>Age category<sup>4</sup></strong></td><td></td><td></td><td></td><td>0,2</td></tr><tr><td>&lt;70</td><td>322 (74%)</td><td>200 (76%)</td><td>122 (71%)</td><td></td></tr><tr><td>≥70</td><td>112 (26%)</td><td>62 (24%)</td><td>50 (29%)</td><td></td></tr><tr><td>Missing</td><td>1</td><td>1</td><td>0</td><td></td></tr><tr><td><strong>Baseline SLD (mm)</strong></td><td>74 (48)</td><td>74 (46)</td><td>75 (51)</td><td>0,9</td></tr><tr><td>Missing</td><td>18</td><td>13</td><td>5</td><td></td></tr><tr><td><strong>Histology<sup>5</sup></strong></td><td></td><td></td><td></td><td>0,042</td></tr><tr><td>Non-squamous</td><td>311 (71%)</td><td>196 (75%)</td><td>115 (67%)</td><td></td></tr><tr><td>Other</td><td>52 (12%)</td><td>33 (13%)</td><td>19 (11%)</td><td></td></tr><tr><td>Squamous</td><td>72 (17%)</td><td>34 (13%)</td><td>38 (22%)</td><td></td></tr><tr><td><strong>Biopsy site</strong></td><td></td><td></td><td></td><td>0,2</td></tr><tr><td>metastasis</td><td>150 (37%)</td><td>95 (39%)</td><td>55 (33%)</td><td></td></tr><tr><td>prim. tumor</td><td>260 (63%)</td><td>148 (61%)</td><td>112 (67%)</td><td></td></tr><tr><td>Missing</td><td>25</td><td>20</td><td>5</td><td></td></tr><tr><td><strong>Treatment line<sup>6</sup></strong></td><td></td><td></td><td></td><td>&lt;0,001</td></tr><tr><td>1</td><td>268 (62%)</td><td>188 (71%)</td><td>80 (47%)</td><td></td></tr><tr><td>2+</td><td>167 (38%)</td><td>75 (29%)</td><td>92 (53%)</td><td></td></tr><tr><td><strong>Liver metastases</strong></td><td></td><td></td><td></td><td>0,11</td></tr><tr><td>FALSE</td><td>369 (85%)</td><td>229 (87%)</td><td>140 (81%)</td><td></td></tr><tr><td>TRUE</td><td>66 (15%)</td><td>34 (13%)</td><td>32 (19%)</td><td></td></tr><tr><td><strong>ECOG<sup>7</sup></strong></td><td></td><td></td><td></td><td>0,006</td></tr><tr><td>0</td><td>174 (40%)</td><td>119 (45%)</td><td>55 (32%)</td><td></td></tr><tr><td>1</td><td>261 (60%)</td><td>144 (55%)</td><td>117 (68%)</td><td></td></tr><tr><td><strong>Gender<sup>8</sup></strong></td><td></td><td></td><td></td><td>0,009</td></tr><tr><td>F</td><td>151 (35%)</td><td>104 (40%)</td><td>47 (27%)</td><td></td></tr><tr><td>M</td><td>284 (65%)</td><td>159 (60%)</td><td>125 (73%)</td><td></td></tr><tr><td><strong>Type of drug<sup>9</sup></strong></td><td></td><td></td><td></td><td>&lt;0,001</td></tr><tr><td>Atezolizumab</td><td>52 (12%)</td><td>24 (9,1%)</td><td>28 (16%)</td><td></td></tr><tr><td>Nivolumab</td><td>51 (12%)</td><td>18 (6,8%)</td><td>33 (19%)</td><td></td></tr><tr><td>Pembrolizumab</td><td>64 (15%)</td><td>33 (13%)</td><td>31 (18%)</td><td></td></tr><tr><td>Pembrolizumab + chemotherapy</td><td>268 (62%)</td><td>188 (71%)</td><td>80 (47%)</td><td></td></tr><tr><td><strong>Smoking history<sup>10</sup></strong></td><td></td><td></td><td></td><td>0,15</td></tr><tr><td>current</td><td>154 (36%)</td><td>101 (40%)</td><td>53 (31%)</td><td></td></tr><tr><td>never</td><td>23 (5,4%)</td><td>12 (4,7%)</td><td>11 (6,4%)</td><td></td></tr><tr><td>previous</td><td>248 (58%)</td><td>140 (55%)</td><td>108 (63%)</td><td></td></tr><tr><td>Missing</td><td>10</td><td>10</td><td>0</td><td></td></tr><tr><td><strong>Mutations<sup>11</sup></strong></td><td></td><td></td><td></td><td></td></tr><tr><td>KRAS</td><td>129 (42%)</td><td>85 (45%)</td><td>44 (38%)</td><td></td></tr><tr><td>TP53</td><td>171 (56%)</td><td>110 (59%)</td><td>61 (52%)</td><td></td></tr><tr><td>STK11</td><td>58 (19%)</td><td>30 (16%)</td><td>28 (24%)</td><td></td></tr><tr><td>KEAP1</td><td>52 (17%)</td><td>34 (18%)</td><td>18 (15%)</td><td></td></tr><tr><td><strong>PD-L1 expression</strong></td><td></td><td></td><td></td><td>&lt;0,001</td></tr><tr><td>&lt;1%</td><td>200 (55%)</td><td>96 (45%)</td><td>104 (69%)</td><td></td></tr><tr><td>[1%, 50%]</td><td>77 (21%)</td><td>53 (25%)</td><td>24 (16%)</td><td></td></tr><tr><td>&gt;50%</td><td>85 (23%)</td><td>63 (30%)</td><td>22 (15%)</td><td></td></tr><tr><td>Missing</td><td>73</td><td>51</td><td>22</td><td></td></tr></tbody></table>

Note: 3 monotherapy patients and 1 combination therapy patient excluded because of treatment duration not long enough to be eligible for primary resistance assessment

<footnote n="1">¹ n (%); Mean (SD)</footnote>

<footnote n="2">² SITC 2023 Guidelines</footnote>

<footnote n="3">³ Pearson's Chi-squared test; Wilcoxon rank sum test; Fisher's exact test</footnote>

<footnote n="4">⁴ Age category. Can either be '<70' or '>=70'</footnote>

<footnote n="5">⁵ Histology. Can either be 'Non-squamous', 'Squamous' or 'Other'. Category 'Other' consists of original entries 'UNK: Unknown', 'UC: Undifferentiated carcinoma', 'SCLC: small cell lung carcinoma', 'NEC: neuro-endocrine carcinoma', 'OTHER', 'NSCLC' and missing values.</footnote>

<footnote n="6">⁶ Treatment line of patient included in the study. First line patient were treated with ICI combination therapy. Patient with at least two line of treatments were treated with ICI monotherapy.</footnote>

<footnote n="7">⁷ Performance status ECOG score. Can either be '0' or '1'</footnote>

<footnote n="8">⁸ Sex. Can either be 'M' or 'F'</footnote>

<footnote n="9">⁹ Type of immunotherapy drug used during treatment.</footnote>

<footnote n="10">¹⁰ Smoking history. Can either be 'current', 'previous' or 'never'</footnote>

<footnote n="11">¹¹ Other hidden mutations are: HER2 (16), CTNNB1 (17), DDR2 (14), FBXW7 (13), SMAD4 (17), ERBB4 (27), BRAF (31), FGFR2 (13), FGFR3 (15), PDGFRA (26), PTEN (27), MET (17), TERT (12), NRAS (7), PIK3CA (28), EGFR (13), KIT (21), POLE (26), RET (22), HRAS (1), IDH1 (1), KIF5B (1), ERBB2 (2), ROS1 (19), TTF1 (1), ANTI_CYTOKERATINE (1), CDKN2A (22)</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 1: Workflow of the analysis</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Figure 2: Univariable biostatistics</figcaption></figure>

<figure><img src="image_4.png" /></figure>
<figure><img src="image_5.png" /><figcaption>Figure 3: Multivariable machine learning predictive performances</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 4: Machine learning signature</figcaption></figure>

<figure><img src="image_7.png" /></figure>
<caption>Figure 5: Digital tool to predict individual risk of resistance</caption>

A. Line 1 PrR

<figure><img src="image_8.png" /></figure>

E. Line 1 non-PrR

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

B. Line 2+ PrR

<figure><img src="image_11.png" /></figure>

F. Line 1 non-PrR

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

C. Line 1 PrR

<figure><img src="image_14.png" /></figure>

G. Line 2+ non-PrR

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

D. Line 2+ PrR

<figure><img src="image_17.png" /></figure>

H. Line 2+ non-PrR

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /></figure>
<figure><img src="image_20.png" /><figcaption>Supplementary Figure 1</figcaption></figure>
Footnotes

(a) Advanced NSCLC patients registered for PD1/L1 ICI-based treatment

(b) Informed consent signed

(c) [According to SITC guidelines (Rivzi, 2023 + Tawbi 2023), Best Response according to RECIST v1.1 with] * a **treatment exposure ≥ 6-8 weeks**

(d) * a **treatment exposure ≥ 6 weeks**

(e) * a **treatment exposure > 6 months**

(f) * a **treatment exposure ≥ 6 months**

(g) NCT03833440

(h) Eligible /Reaching biopsy time, N=349
Supplementary Figure 2: Progression Free Survival (PFS)

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /></figure>
Supplementary Figure 2: Progression Free Survival (PFS)

Kaplan-Meier plots of progression-free survival.

A. All patients

B. First-line patients (combination therapy)

C. Subsequent-line patients (monotherapy)
<table><caption>Supplementary Figure 3</caption><thead><tr><th>Signif</th><th>Marker</th><th>AUC</th><th>OR PrR</th><th>HR PFS</th><th>C-Index PFS</th><th>HR OS</th><th>C-Index OS</th></tr></thead><tbody><tr><td rowspan="10">PrR + PFS + OS</td><td>Treg Stroma</td><td>0.64</td><td>0.572 (0.444 - 0.737)</td><td>0.771 (0.679 - 0.877)</td><td>0.59</td><td>0.792 (0.689 - 0.909)</td><td>0.58</td></tr><tr><td>C reactive protein</td><td>0.63</td><td>1.89 (1.45 - 2.47)</td><td>1.52 (1.31 - 1.76)</td><td>0.61</td><td>1.85 (1.54 - 2.21)</td><td>0.65</td></tr><tr><td>Chloride</td><td>0.63</td><td>0.616 (0.478 - 0.794)</td><td>0.685 (0.593 - 0.792)</td><td>0.61</td><td>0.604 (0.512 - 0.712)</td><td>0.63</td></tr><tr><td>Albumin</td><td>0.60</td><td>0.616 (0.481 - 0.789)</td><td>0.727 (0.634 - 0.833)</td><td>0.58</td><td>0.644 (0.555 - 0.747)</td><td>0.62</td></tr><tr><td>PD-L1</td><td>0.62</td><td>0.672 (0.533 - 0.848)</td><td>0.661 (0.566 - 0.772)</td><td>0.60</td><td>0.677 (0.562 - 0.815)</td><td>0.58</td></tr><tr><td>Red cells</td><td>0.61</td><td>0.628 (0.496 - 0.796)</td><td>0.673 (0.586 - 0.774)</td><td>0.61</td><td>0.61 (0.519 - 0.716)</td><td>0.61</td></tr><tr><td>VCAM</td><td>0.64</td><td>2.15 (1.36 - 3.4)</td><td>1.34 (1.16 - 1.56)</td><td>0.60</td><td>1.53 (1.29 - 1.83)</td><td>0.62</td></tr><tr><td>ICAM</td><td>0.63</td><td>1.97 (1.32 - 2.93)</td><td>1.46 (1.19 - 1.8)</td><td>0.59</td><td>1.48 (1.17 - 1.86)</td><td>0.60</td></tr><tr><td>ECOG</td><td>0.57</td><td>1.99 (1.27 - 3.12)</td><td>1.68 (1.29 - 2.19)</td><td>0.57</td><td>1.78 (1.31 - 2.42)</td><td>0.56</td></tr><tr><td>ASAT</td><td>0.57</td><td>1.47 (1.17 - 1.85)</td><td>1.25 (1.1 - 1.42)</td><td>0.56</td><td>1.23 (1.07 - 1.4)</td><td>0.55</td></tr><tr><td>Platelets</td><td>0.56</td><td>1.45 (1.14 - 1.84)</td><td>1.25 (1.1 - 1.43)</td><td>0.54</td><td>1.37 (1.18 - 1.58)</td><td>0.57</td></tr><tr><td rowspan="4">PrR + OS</td><td>Sodium</td><td>0.57</td><td>0.689 (0.545 - 0.87)</td><td>0.773 (0.675 - 0.884)</td><td>0.57</td><td>0.681 (0.585 - 0.794)</td><td>0.59</td></tr><tr><td>Alkaline phosphatase</td><td>0.55</td><td>1.56 (1.21 - 2)</td><td>1.46 (1.27 - 1.68)</td><td>0.56</td><td>1.41 (1.26 - 1.59)</td><td>0.60</td></tr><tr><td>NLR</td><td>0.57</td><td>1.46 (1.14 - 1.87)</td><td></td><td></td><td>1.38 (1.16 - 1.64)</td><td>0.58</td></tr><tr><td>PD-L1 expression 1%</td><td></td><td></td><td>1.53 (1.12 - 2.09)</td><td>0.58</td><td>1.68 (1.16 - 2.44)</td><td>0.57</td></tr><tr><td rowspan="8">PFS + OS</td><td>Transitional B cells</td><td></td><td></td><td>0.794 (0.694 - 0.909)</td><td>0.58</td><td>0.726 (0.626 - 0.842)</td><td>0.63</td></tr><tr><td>Activated T cells (%)</td><td></td><td></td><td>1.23 (1.08 - 1.39)</td><td>0.57</td><td>1.25 (1.08 - 1.44)</td><td>0.56</td></tr><tr><td>CD8 ILT2+ T cells (%)</td><td></td><td></td><td>1.21 (1.05 - 1.38)</td><td>0.57</td><td>1.24 (1.07 - 1.45)</td><td>0.55</td></tr><tr><td>IgA+ B cells (%)</td><td></td><td></td><td>1.24 (1.07 - 1.43)</td><td>0.53</td><td>1.36 (1.15 - 1.61)</td><td>0.54</td></tr><tr><td>MLR</td><td></td><td></td><td>1.28 (1.13 - 1.46)</td><td>0.57</td><td>1.41 (1.22 - 1.64)</td><td>0.59</td></tr><tr><td>MDSC M MDSC Stroma Tumor</td><td></td><td></td><td>0.831 (0.732 - 0.944)</td><td>0.56</td><td>0.828 (0.718 - 0.955)</td><td>0.57</td></tr><tr><td>Liver metastases</td><td></td><td></td><td>1.64 (1.18 - 2.28)</td><td>0.53</td><td>1.86 (1.3 - 2.66)</td><td>0.55</td></tr><tr><td>Dendritic cells</td><td></td><td></td><td>1.17 (1.04 - 1.31)</td><td>0.50</td><td>1.21 (1.08 - 1.35)</td><td>0.51</td></tr><tr><td rowspan="4">PFS</td><td>Ferritin</td><td></td><td></td><td>1.49 (1.14 - 1.96)</td><td>0.59</td><td>1.69 (1.2 - 2.39)</td><td>0.61</td></tr><tr><td>Marginal zone-like B cells</td><td></td><td></td><td>0.845 (0.751 - 0.951)</td><td>0.57</td><td>0.819 (0.726 - 0.924)</td><td>0.58</td></tr><tr><td>Treg cells</td><td></td><td></td><td>0.802 (0.691 - 0.93)</td><td>0.52</td><td>0.789 (0.694 - 0.897)</td><td>0.57</td></tr><tr><td>Heart frequency (bpm)</td><td></td><td></td><td>1.21 (1.05 - 1.38)</td><td>0.54</td><td>1.25 (1.06 - 1.46)</td><td>0.53</td></tr><tr><td rowspan="8">OS</td><td>PD-L1+ cell dens. in Tumor</td><td></td><td></td><td>0.827 (0.731 - 0.935)</td><td>0.58</td><td></td><td></td></tr><tr><td>PD1+ cells</td><td></td><td></td><td>0.72 (0.602 - 0.862)</td><td>0.61</td><td></td><td></td></tr><tr><td>CD3+CD4+ stroma</td><td></td><td></td><td>0.817 (0.728 - 0.917)</td><td>0.59</td><td></td><td></td></tr><tr><td>PD-L1+ non classical monocytes</td><td></td><td></td><td>1.28 (1.14 - 1.43)</td><td>0.52</td><td></td><td></td></tr><tr><td>Presence of mutations on ANTI-CYTOKERATINE</td><td></td><td></td><td>38.3 (4.4 - 333)</td><td>0.51</td><td></td><td></td></tr><tr><td>Presence of mutations on CTNNB1</td><td></td><td></td><td>2.49 (1.44 - 4.31)</td><td>0.53</td><td></td><td></td></tr><tr><td>Number of metastases at baseline</td><td></td><td></td><td>1.2 (1.06 - 1.36)</td><td>0.53</td><td></td><td></td></tr><tr><td>Prothrombin rate</td><td></td><td></td><td></td><td></td><td>0.777 (0.682 - 0.885)</td><td>0.60</td></tr><tr><td>Presence of mutations on STK11</td><td></td><td></td><td></td><td></td><td>1.68 (1.14 - 2.46)</td><td>0.56</td></tr><tr><td>ISIC Three Groups</td><td></td><td></td><td></td><td></td><td>2.08 (1.27 - 3.4)</td><td>0.56</td></tr><tr><td>ISIC CD8free d.20 tumor</td><td></td><td></td><td></td><td></td><td>0.793 (0.689 - 0.914)</td><td>0.58</td></tr><tr><td>CD69+ HLA-DR+ B cells</td><td></td><td></td><td></td><td></td><td>0.784 (0.671 - 0.916)</td><td>0.56</td></tr><tr><td>LDH</td><td></td><td></td><td></td><td></td><td>1.36 (1.15 - 1.61)</td><td>0.60</td></tr><tr><td>M MDSC</td><td></td><td></td><td></td><td></td><td>1.23 (1.08 - 1.4)</td><td>0.58</td></tr><tr><td>Histology</td><td></td><td></td><td></td><td></td><td>1.66 (1.17 - 2.36)</td><td>0.54</td></tr><tr><td>CD4 CD39+ T cells</td><td></td><td></td><td></td><td></td><td>0.812 (0.695 - 0.948)</td><td>0.56</td></tr><tr><td>Serum sCD25</td><td></td><td></td><td></td><td></td><td>1.23 (1.07 - 1.41)</td><td>0.56</td></tr><tr><td>Basic Mass Index</td><td></td><td></td><td></td><td></td><td>0.792 (0.671 - 0.935)</td><td>0.55</td></tr><tr><td>Serum TNFa</td><td></td><td></td><td></td><td></td><td>0.791 (0.675 - 0.927)</td><td>0.53</td></tr></tbody></table>
<figure><img src="image_24.png" /><figcaption>Supplementary Figure 4</figcaption></figure>