Distinct immune cell dynamics associated with immune-related adverse events during combined chemoradiation and immune checkpoint inhibitor therapy

Domain: chemistry
Published: 2025-12-22
Editorial outcome: major_revision

Abstract
Combining chemoradiotherapy with immunotherapy increases the risk of immune-related adverse events (irAEs), but the underlying mechanisms remain poorly understood. To address this, we conducted a longitudinal single-cell multi-omics analysis of patients with locally advanced cervical cancer. Here we show that the proportions of CD4+ and CD8+ terminally differentiated effector memory or effector T cells are elevated in patients with irAEs. Chemoradiotherapy reduces B cell clonality while increasing the abundance and somatic hypermutation frequencies of IgA+ and IgG+ B cells in irAE patients. In the myeloid compartment, combined treatment expands specific monocyte subclusters associated with irAEs. Spatial transcriptomics and immunofluorescence analyses further reveal that these irAE-associated immune cells aggregate within the tumor microenvironment. Finally, we develop the predictive models for irAEs and integrate them, along with all datasets, into a user-friendly data portal. Our findings suggest that chemoradiotherapy and immunotherapy exert distinct effects on different immune cells, contributing to irAE development.

Introduction
Cervical cancer is the third most common cause of death among women globally1, with its incidence continuing to rise annually in China2. Although early-stage cervical cancer is generally curable, approximately 40–50% of patients are diagnosed with locally advanced cervical cancer (LACC, International Federation of Gynecology and Obstetrics stage IB3-IVA), and this contributes significantly to the high mortality associated with this disease3,4. The current standard treatment for LACC is concurrent chemoradiotherapy. However, 30–50% of these patients experience disease recurrence or progression within five years5,6. Combining chemoradiotherapy with immunotherapy is being explored as a potential strategy to enhance treatment efficacy and improve the prognosis for patients with LACC7. However, this combination also brings an increased risk of immune-related adverse events (irAEs), as the immune responses triggered by radiotherapy may synergize with those from immunotherapy8,9. In addition to divergent treatment efficacies observed in recent clinical trials, the incidence of grade 3 or higher treatment-related toxicities remains consistently high5,6,10–12.

IrAEs are caused by excessive activation of the immune system against host tissues, potentially leading to serious consequences, such as discontinuation of treatment, permanent tissue damage and even fatal outcomes13,14. These reactions, which vary from mild to life-threatening, can significantly diminish the clinical benefits to the patients15. Nevertheless, a few studies have identified molecular characteristics associated with treatment efficacy rather than irAEs. For example, enhanced stromal immune infiltration10, higher tumor aneuploidy scores16 and the presence of circulating antigen-experienced CD103+ T cells17 have all been linked to improved patient outcomes. Notably, immune responses specifically related to irAEs in this context remain largely unexplored. Understanding the mechanisms driving irAEs and identifying relevant biomarkers could improve the prediction and management of these adverse events, ultimately enhancing treatment safety and clinical outcomes for patients undergoing combination immunotherapy.

To address existing knowledge gaps, we conducted a clinical trial (NCT06378840, an observational standard of care treatment study) aimed to evaluate the role of immune cells in treatment-naive patients with LACC receiving sequential chemoradiotherapy and ICI therapy. We performed single-cell transcriptome and V(D)J sequencing on matched tumor and peripheral blood samples, which were collected at four distinct treatment stages. This approach allowed us to investigate how different treatments affected both tumor behavior and systemic immune responses, with a specific focus on identifying irAE-specific changes triggered by various treatment approaches. By utilizing these longitudinal data, our deep learning models was able to identify representative irAE characteristics of immune cells at various treatment stages, and this was further validated in ten independent irAE datasets. This study provides a comprehensive single-cell immune profile related to irAEs and establishes reliable predictive models for irAEs. These models may help identify patients with higher risks of developing irAEs and could potentially assist in optimizing patient selection for those undergoing combined chemoradiation and ICI therapies.

Results
Longitudinal single-cell multi-omics atlas of LACC tumor and blood samples
We prospectively enrolled thirty patients with treatment-naive LACC and collected 96 tumor biopsies and peripheral blood samples. Due to the prospective nature of our study and the need to maintain high cell integrity and RNA quality, we exclusively used fresh samples and did not include any frozen samples. Among the patients, twenty received sequential chemoradiotherapy combined with anti-programmed cell death-1 (PD-1) antibodies (radio+ICI group), while the remaining ten underwent chemoradiotherapy followed by chemotherapy alone (radio group) (Fig. 1a, b, Supplementary Table 1a). All patients were closely monitored for irAE development throughout their treatment, with a median follow-up time of 524 days and a median time to irAE of 139 days (Supplementary Table 1b). The majority of patients experienced at least one irAEs, ranging from mild (grade 1) to severe (grade 3), and these affected various organ systems (Supplementary Table 1b, Methods). Blood and tumor samples were collected at key time points: pre-chemoradiotherapy (pre-radio), during chemoradiotherapy (on-radio), pre- and during ICI treatment (pre-imm and on-imm) for radio+ICI group and pre- and during chemotherapy (pre-chemo and on-chemo) for radio group (Fig. 1a). Specifically, pre-radio samples were collected before initiation of therapy. On-radio samples were collected after the 10th fraction of external beam irradiation (1.8 Gy per fraction). Pre-imm and pre-chemo samples were obtained after completing chemoradiotherapy but before administering ICI or chemotherapy. On-imm samples were collected before the 4th cycle of ICI, and on-chemo samples were taken before the 3rd cycle of chemotherapy (Fig. 1a). This sample collection strategy allowed us to track local and systemic immunological dynamics throughout the treatment process.

We employed 10× Genomics Chromium single-cell RNA sequencing (scRNA-seq) in conjunction with single-cell V(D)J sequencing (scV(D)J-seq) of T cell receptor (TCR) and B cell receptor (BCR) clonotypes to analyze the molecular characteristics and alterations of the tumor microenvironment (TME) and systemic immune system in patients with LACC during treatment (Fig. 1b, Supplementary Table 1a). After quality control and removal of doublets, the scRNA-seq assay yielded 708,924 cells, representing 13 major cell types and 48 fine-grained cell subclusters based on the canonical marker gene expression (Fig. 1c, Supplementary Fig. 1a). These cell types were shared across patients, tissue types, sample collection time points and treatment modalities (Supplementary Fig. 1b–e). We identified the tumor cells from epithelial cells by using the inferred copy number variation (CNV), which aligned with the data obtained from scRNA-seq and whole exome sequencing data from another study18, with the cells generally being clustered by patient (Supplementary Fig. 1b, f).

Cellular dynamics in response to chemoradiotherapy and ICI treatment in relation to irAEs
Next, we analyzed the cellular abundance of major cell types in relation to treatments and clinical endpoints, defined by the occurrence of irAEs and disease relapse. Chemoradiotherapy significantly reduced tumor cell proportions in tumor samples and increased macrophages in both tumor and peripheral blood samples, with monocytes showing a notable increase exclusively in peripheral blood samples (Supplementary Fig. 2a, b), consistent with well-documented effects of radiotherapy19,20. Following ICI treatment, T cells, macrophages and monocytes showed slight decreases in blood samples (Supplementary Fig. 2b). Comparisons between irAE and non-irAE patients revealed elevated proportions of B cells, T cells and macrophages in blood samples collected from irAE patients prior to ICI administration. After ICI treatment, T cells increased, while macrophages decreased in irAE patients (Fig. 1d, Supplementary Fig. 3d, e). No significant differences in immune cell proportions were observed in tumor samples at pre-radio and on-radio time points between the irAE and non-irAE patients (Supplementary Fig. 3a–b). Collectively, these findings indicate that chemoradiotherapy reduces tumor cells and increases macrophages across tumor and peripheral blood samples. Notably, irAE-related immune cell changes were observed in the peripheral blood samples collected before and after ICI treatment, suggesting that systemic immune responses appear to be more closely linked to irAEs. However, pre-radio tumor samples from patients who remained relapse-free two years after receiving first therapy (Methods) showed increased proportions of several immune cells, including B cells, T cells and macrophages. Meanwhile, blood samples from non-relapsed patients collected prior to ICI treatment displayed elevated levels of B cells and T cells, with macrophages decreased before and after ICI treatment (Fig. 1e, Supplementary Fig. 3f–i). These results indicate that treatment response-related alterations were predominantly identified in baseline tumor samples and immune responses within the TME may be more strongly associated with treatment outcomes.

Increased CD8.Temra and CD4.Temra cells during combined treatment are associated with irAEs
Our aforementioned analyses revealed an increased proportion of T cells in the blood and tumor associated with irAEs and disease relapse, respectively. Next, we compared the abundances and functions of T cell subclusters between patients in different irAE and disease relapse subgroups. We used canonical immune markers21 and curated gene signatures22 to identify T cell subclusters and these included eight CD4+, eight CD8+, two natural killer T (NKT) as well as one each of γδ T and proliferative T subclusters (Fig. 2a–d, Supplementary Figs. 4a, 5a–e). While most T cell subclusters were shared across samples from different patients as well as during treatments and at different time points, some of these, such as CD4+ terminally differentiated effector memory or effector T cells (CD4.Temra), CD8+ terminally differentiated effector memory or effector T cells (CD8.Temra) and CD4+ adenylosuccinatelyase-positive naive T cells (CD4.ADSL+ Tn), showed a strong distribution preference in peripheral blood (Fig. 2a, Supplementary Fig. 4b). For CD4+ T cells in patients with irAEs, the proportions of CD4.Terma, with cytotoxicity phenotype, and CD4.ADSL+ Tn, with a highly expressed oxidative phosphorylation signature, were significantly reduced in baseline pre-radio blood samples but increased in pre-imm blood samples after chemoradiotherapy. Consistent with the reported role of memory CD4+ T cells in irAEs23, CD4.Temra increased following ICI treatment in both irAE groups, with higher abundances observed in patients with irAEs (Fig. 2e, f, Supplementary Fig. 5e).

For CD8+ T cell subclusters, we observed significantly increased proportions of CD8.Temra, with highly expressed cytotoxicity signature, and CD8+ granzyme K–positive effector memory T cells (CD8.GZMK+ Tem), with highly expressed senescence and stress response signatures, in both the pre-radio (baseline) and pre-imm (after chemoradiotherapy but before ICI treatment) blood samples of irAE patients. However, these two subclusters showed opposite alterations in on-imm (after chemoradiotherapy and during ICI treatment) blood samples from irAE patients (Fig. 2e, f, Supplementary Fig. 5d). Interestingly, two NK-like CD8+ T cell subclusters, with relatively higher expression of cytotoxicity signature (Supplementary Fig. 5d), and one NKT subcluster, CD8+ killer cell immunoglobulin-like receptor-positive tyrosine kinase-positive NK-like T cell (CD8.KIR+ TXK+ NK-like), CD8+ killer cell immunoglobulin-like receptor-positive eomesodermin-negative NK-like T cell (CD8.KIR+ EOMES-NK-like) and fibroblast growth factor binding protein 2-positive NK T (NK T.FGFBP2+), decreased in both pre-imm and on-imm samples (Fig. 2e). Our findings indicated that during chemoradiotherapy and ICI treatment, the proportions of CD4.Temra and CD8.Temra, with cytotoxicity and memory phenotypes, were elevated in patients experiencing irAEs, indicating the involvement of both CD4+ and CD8+ T cells in irAE development.

Although the overall T cell proportions were significantly higher in the pre-radio tumor sample of non-relapse patients, we only observed proportion increments in three CD4+ T cell subpopulations. These included CD4+ capping actin protein, gelsolin like-positive memory T cells (CD4.CAPG+ Tm), CD4+ tumor necrosis factor receptor superfamily member 9-positive regulatory T cells (CD4.TNFRSF9+ Treg) and CD4.Terma (Supplementary Fig. 5f), which exhibited highly expressed cytotoxicity or exhaustion signatures (Supplementary Fig. 5e).

Differential TCR diversity and clonal expansion in CD4+ and CD8+ T cells of irAE patients following ICI treatment
Using TCR sequences (clonotypes) as molecular barcodes enables mapping of T cell clonal expansions and diversity resulting from chemoradiotherapy and ICI administration. We obtained 161,711 cells with both α and β TCR chains from the CD4+ and CD8+ subsets, including 79,811 with unique and 81,900 with repeated TCRs (Supplementary Fig. 6a–d). The repeated TCRs are indicative of clonal expansion. TCR diversity, assessed by the Shannon indexes, did not show a significant correlation with disease relapse for either CD8+ or CD4+ T cells (Supplementary Fig. 6e), suggesting limited disease relapse-related differences in TCR within current dataset. Therefore, subsequent analysis focused on irAEs. Notably, CD4+ T cells from irAE patients displayed significantly higher TCR diversity in on-imm blood samples compared to non-irAE patients, with no significant differences observed at the other time points (Fig. 3a). Conversely, the TCR diversity of CD8+ T cells shows no significant association with irAEs (Supplementary Fig. 6f). These findings suggest that irAE patients have a more polyclonal repertoire of CD4+ T cells but not CD8+ T cells. Clonal expansions were observed in both CD4+ and CD8+ T cells across blood and tumor samples, with divergent percentages of overlapping clones (0.47 to 9.02% and 4.42 to 41.48% for CD4+ and CD8+ T cells, respectively; Fig. 3b, Supplementary Fig. 6g). In alignment with the decreased TCR diversity, overlapping CD4+ T cell clonotypes between pre-imm and on-imm blood samples, as well as between pre-imm and pre-radio blood samples, were significantly lower in irAE patients compared to non-irAE patients. In contrast, no significant differences were observed in overlapping CD8+ T cells clonotypes (Fig. 3b, Supplementary Fig. 6g).

Given the high fractions of expanded CD8+ T cells, we next investigated whether TCRs from clonally expanded CD8+ T cells were shared between subclusters that were increased in irAE patients. In irAE patients, CD8.Temra and CD8.GZMK+ Tem exhibited significantly higher percentages of clonally expanded TCRs shared with CD8.Temra after ICI treatment (Fig. 3c). However, no significant changes were observed in CD4+ T cells (Supplementary Fig. 6h). We next asked whether the expanded TCRs recognized specific epitopes. As HPV abundance was negatively associated with the proportions of most CD8+ T subclusters (Fig. 3d, Supplementary Fig. 6i), we explored whether CD8+ T cells harbored TCRs recognizing HPV epitopes. Interestingly, six irAE patients had HPV-specific TCRs in samples collected at all time points, whereas only three non-irAE patients had HPV-specific TCRs in samples collected at baseline (Fig. 3e). These results suggest that expanded HPV-recognizing T cells may play a role in irAE development. In summary, our findings indicate that CD4+ T cells displayed greater TCR diversity and fewer overlapping clonotypes across time points and tissues obtained from irAE patients. Meanwhile, CD8+ T cells exhibited a higher proportion of TCRs shared with irAE-associated subclusters after ICI treatment, suggesting distinct roles of CD4+ and CD8+ T cells in irAEs.

Lower BCR clonality after chemoradiotherapy in irAE patients
Previous studies have reported associations between B cells and irAEs, such as clonality24 and abundance of specific B cells25. Based on this, we investigated the abundance and immunoglobulin (Ig) repertoires of B-lineage cells in patients with and without irAEs. Unsupervised clustering identified five distinct B cell subclusters: two naive B cell subclusters (tubulin beta 1 class VI-positive naive B cell, Bn-TUBB1, and T-cell leukemia/lymphoma protein 1A-positive naive B cell, Bn-TCL1A), one memory B cell subcluster, one transitional B subcluster and one plasma cell (PC) subcluster (Fig. 4a, Supplementary Fig. 7a). Unlike T cells, B cell subclusters exhibited few significant frequency changes across patient groups. Proportions of Bn-TCL1A significantly increased following ICI treatment, especially in irAE patients (Fig. 4b, Supplementary Fig. 7b, c). The percentages of the top BCR clonotypes at on-radio, pre-imm and on-imm time points increased in non-irAE patient but not in irAE patients, indicating a more even distribution of clonal size in irAE patients following chemoradiotherapy and ICI treatment (Fig. 4c). Next, we analyzed Ig heavy-chain isotypes based on the Ig constant region sequences. As expected, naive B and transitional B cell subclusters were dominated by IgM isotypes, memory B cell exhibited a variety of isotypes and the PCs were highly enriched in IgG and IgA isotypes (Supplementary Fig. 8a, b). Comparisons between irAE and non-irAE patients revealed significant increases in the proportion of IgA+ and IgG+ B cells in both blood and tumor samples following chemoradiotherapy, accompanied by a decrease in IgM+ B cells (Fig. 4d). To further investigate the functional status of B cells, we analyzed somatic hypermutation (SHM), where high SHM levels indicate strong antigen-specific affinity26. As expected, most naive B cells exhibited no SHM, memory and transitional B cells displayed moderate SHM frequencies, and PCs had higher SHM frequencies (Fig. 4e). Consistent with the observed isotype dynamics, SHM frequencies of IgA+ and IgG+ B cells increased in both blood and tumor samples from irAE patients following chemoradiotherapy (Fig. 4f). Overall, these findings indicate that chemoradiotherapy reduced B cell clonality while increasing the abundance and SHM frequencies of IgA+ and IgG+ B cells in blood and tumor samples of patients with irAEs. This dynamic may play a role in the formation of autoantibodies, ultimately leading to irAEs.

Analyses of disease relapse status showed that memory B cells were more abundant in pre-radio tumor samples from patient without disease relapse (Supplementary Fig. 7d). However, the distribution of BCR clonotypes were similar between patients with and without relapse (Supplementary Fig. 8c). IgA+ and IgM+ B cells exhibited inverse changes in tumor samples from non-relapsed patients after chemoradiotherapy (Supplementary Fig. 8d). SHM frequencies remained largely unchanged in relapsed and non-relapsed patients (Supplementary Fig. 8e). Collectively, although memory B cell increased in non-relapsed patients, clonotypes and SHM did not change significantly.

Divergent myeloid cell reprogramming induced by a combination of chemoradiation and ICI therapies
Building on our previous study suggesting the distinct roles of macrophages in the TME related to irAEs and treatment efficacies in lung squamous cell carcinoma15, we further investigated myeloid cells. We identified five major myeloid cell populations, further categorized into 14 unique subclusters based on published signatures and cell markers27 (Fig. 5a–d, Supplementary Figs. 9, 10). The frequencies of myeloid cell subclusters differed significantly between irAE and non-irAE groups across various time points. Notably, two C-X3-C motif chemokine receptor 1 (CX3CR1)-positive macrophage subclusters displayed contrasting trends at the pre-imm time point. Tumor necrosis factor receptor superfamily member 8-positive macrophage (Macro.TNFRSF8), characterized by a high inflammatory response signature, significantly increased in irAE patients, whereas G protein-coupled receptor kinase 2-positive macrophages (Macro.GRK2) decreased (Fig. 5e, Supplementary Fig. 9d). Among monocyte subclusters, folate receptor 3-positive (Mono.FOLR3) monocyte, which expressed the highest inflammatory response signature, were enriched in pre-imm and on-imm blood samples from irAE patients (Fig. 5e). Conversely, four monocyte subclusters with lower expression of inflammatory response, C-C motif chemokine ligand 5-positve (Mono.CCL5), lymphocyte antigen 96-positive (Mono.LY96), serpin family B member 1-positive (Mono.SERPINB1) and the major histocompatibility complex class II-positive (Mono.MHCII) monocyte, were reduced in irAE patients both before and after ICI treatment (Fig. 5e). When analyzing gene signature expression related to myeloid cell functions, most macrophage and monocyte subclusters in irAE patients’ pre-imm blood samples exhibited significantly higher levels of inflammation- and antigen presentation-related pathways (Fig. 5f). In contrast, at the on-imm time point, only subclusters with increased frequencies in irAE patients maintained elevated inflammation- and angiogenesis-related signature expression. However, phagocytosis-related signatures were consistently higher in most myeloid subclusters in irAE patients (Fig. 5f). Other gene signatures expressed varied across samples collected at pre-radio and on-radio time points (Supplementary Fig. 10h). In summary, chemoradiotherapy combined with ICI treatment increased the frequencies of monocyte subclusters with heightened inflammatory responses and enhanced the phagocytic capabilities of myeloid cells in patients with irAEs. These findings suggest that the inflammatory status of monocytes may contribute to irAE development.

Two macrophage subclusters with higher phagocytosis-related signature expression, tryptase alpha/beta 1-positive (Macro.TPSAB1) and triggering receptor expressed on myeloid cells 2-positive (Macro.TREM2) macrophage, were enriched in baseline tumor samples from non-relapsed patients, but interleukin 1 beta-positive macrophage (Macro.IL1B), with higher expression of inflammatory response, was reduced (Supplementary Figs. 9d, 10g). Mono.CCL15 and Mono.LY96 decreased in pre-imm blood samples of non-relapsed patients (Supplementary Fig. 10g). Myeloid cell function-related analysis found that in non-relapsed patients, angiogenesis-related signatures were elevated in most macrophage and monocyte subclusters both before and after chemoradiotherapy and ICI treatment (Supplementary Fig. 10i). Altogether, these results suggested that angiogenic potential of myeloid cells may positively influence treatment efficacy.

Colocalization of cells associated with irAEs within the TME in response to chemoradiotherapy
We performed Spearman correlation analysis on the frequencies of cell subclusters in pre-radio and on-radio tumors in order to explore the relationships between cellular dynamics related to irAEs. In the on-radio tumor samples, we found strong correlations among cell subclusters that showed significant changes in samples collected after the on-radio time point, such as in the CD4.Temra, CD8.Temra and CD4.ADSL+ Tn (Supplementary Fig. 11a). However, the correlation coefficients were similar in the pre-radio tumor samples (Supplementary Fig. 11b). In order to further investigate the spatial organization of irAE-associated cells within the TME of irAE patients, we performed spatial transcriptomic analysis using the 10× Genomics Xenium platform, which profiled 380 genes. Subcellular resolution of the Xenium data28 allowed us to capture and visualize the spatial localization of different cell subclusters within the TME (Supplementary Fig. 12). We observed a significant decrease in malignant cells and an increase in fibroblasts in on-radio tumors when compared to the pre-radio tumors, and this aligned with the results from the typical radiation-induced changes observed in the tumors29 (Fig. 6). Notably, the immune cells associated with irAEs, including B cell, CD4+ T cells and CD8+ T cells, were aggregated in the on-radio tumors of irAE patients, whereas this aggregation was absent in the on-radio tumors obtained from non-irAE patients (Fig. 6a, c). To corroborate this observation, we performed the multiplex immunofluorescence (mIF) staining using markers for total T cells (CD3), CD4+ T cells (CD4), CD8+ T cells (CD8), B cells (CD20), macrophages (CD68) and malignant cells (pan-CK) (Fig. 6e). We then analyzed the spatial distribution of immune cells and measured the intercellular distances between B cells, CD4+ T cells and CD8+ T cells using HALO software (Fig. 6f). The results indicated an increased proximity between immune cells in irAE patients compared to non-irAE patients. Furthermore, we quantified the number of CD20+ in close proximity to CD4+, CD8+ and CD3+ cells (distance <= 50 μm) and observed significantly more such cell proximities in tumor samples from irAE patients (Fig. 6g), suggesting increased immune cell aggregations and interactions in the TME. This spatial organization of key immune cells may trigger subsequent immune responses, potentially contributing to the development of irAEs.

The predictive model based on the molecular characteristics of immune cells during ICI treatment showed superior performance in predicting irAEs
Given that the majority of irAE-related alterations were observed in peripheral blood samples and that different cell types exhibited divergent molecular features, we aimed to develop predictive models for irAEs using the scRNA-seq data from pre-radio, pre-imm and on-imm blood samples. Initially, we employed a method that based on differentially expressed genes (DEG) in cell subclusters enriched in irAE patients (DEG-based irAE signatures, Methods). However, we found that none of the signatures consistently showed increased levels across the validation cohort, regardless of whether scRNA-seq or bulk RNA sequencing (RNA-seq) datasets were used (Supplementary Fig. 13). Cumulative incidence analysis suggested that low irAE scores indicated higher irAE risk or failed to achieve statistical significance (Supplementary Fig. 13f, g). These results suggest that DEG-based irAE signatures may have limitations in capturing the full molecular complexity of irAEs. Therefore, we explored deep learning strategies to enhance irAE prediction. We developed models based on the scRNA-seq data obtained from pre-radio, pre-imm, and on-imm blood samples, referred to as models 1, 2, and 3, respectively. Immune cell subclusters at different time points were classified into irAE-related, significantly enriched in irAE patients, and non-irAE-related groups, significantly enriched in non-irAE patients. Genes were ranked for both groups based on their expressions in individual cells, and ‘cell sentences’ were formulated using gene names and their ranks. These ‘cell sentences’ were then input into the bidirectional encoder representation from transformer (BERT) language model, forming the foundation of our irAE prediction models (Methods). We performed subsampling analyses to assess the robustness of these models. We randomly partitioned these labeled cells into training sets consisting of 70%, 80% and 90% of the data, with corresponding testing sets of 30%, 20% and 10%, respectively. Model performance was evaluated using area under the receiver operating characteristic curve (AUC), Precision-Recall-AUC (PRAUC), accuracy (ACC) and F1 score, based on the testing data. Across all four metrics, the models demonstrated stable performance (Supplementary Fig. 14a), indicating robustness.

To ensure these models accurately captured representative features of irAEs at various treatment stages, we validated them using samples (Supplementary Table 2) collected at time points comparable to those used for corresponding model development. Specifically, we used blood samples collected prior to ICI (pre-ICI blood) to validate models 1 and 2, blood samples collected during ICI treatment (on-ICI blood) to validate model 3 and their AE-affected tissue, where almost all irAE-related immune responses occurred, to validate all models. In pre-ICI blood from two independent datasets30,31, the model 2, model 1, consistently predicted higher irAE scores in patients who developed irAEs (Fig. 7a, Supplementary Fig. 14b). In three additional datasets31–33 with on-ICI blood samples, the predicted irAE scores from the model 3 were also significantly higher in irAE patients (Fig. 7a, Supplementary Fig. 14b). Interestingly, in four datasets collected tissues affected by irAEs33–36, such as colitis and myocarditis, both the model 2 and 3 consistently predicted higher irAE scores in irAE-affected tissues compared to control samples (Fig. 7b, Supplementary Fig. 14c). These findings demonstrate that the deep learning algorithm effectively captured irAE-related molecular features across different treatment time point.

To broaden the applicability of irAE models and identify the most effective model, we applied them to bulk RNA-seq data obtained from blood samples collected before ICI administration. The time points before treatment are critical for irAE prediction, as these stages provide valuable opportunities to implement interventions or alternative treatment designs to benefit patients. Our analysis included an in-house cohort (n=14, Supplementary Table 3) with pre-imm blood samples from patients treated with combined radiotherapy and ICI, as well as two publicly available datasets30,37 with larger sample sizes (n=53 and n=243, respectively). In these three datasets, the model 2 and 3 predicted significantly higher irAE scores in pre-ICI blood samples from irAE patients, with the model 3 showing the most pronounced differences (Fig. 7c, d). Cumulative incidence analyses of the time to any irAE onset37 or severe irAE onset30, based on the follow-up duration after ICI administration, revealed that patients with higher predicted irAE scores from the model 3 had a significantly shorter time to irAE onset, with hazard ratios of 1.50 and 2.90, respectively, whereas the model 2 did not show the same predictive potential (Fig. 7e, f). These results indicate the predictive potential of the model 3 in forecasting patient susceptibility to irAEs.

Interpreting model predictions via cell subcluster analysis
To better comprehend these models, we performed ablation analyses to quantify the contributions of distinct cell subclusters to model performance. We ablated a single cell subcluster by excluding its data from the input and training the model on the remaining data. The ablated model was then applied to the full dataset to generate predictions, from which we computed AUC, PRAUC, ACC, and F1 score. Consistent with previous observations, in Model 3, ablations of CD8.Temra, Mono.FOLR3 and NK T.FGFBP2 resulted in the top declines across four metrics (Supplementary Fig. 14d–g). In Model 2, removal of Macro.TNFRSF8, CD8.Temra and NK T.FGFBP2 caused the top reductions across these metrics (Supplementary Fig. 14d–g). Next, we exported the top 100 genes (Supplementary Table 4, determined by keyword extraction approach for BERT model38) from model 2 and model 3 and found that these genes were enriched in pathways associated with activated immune responses, such as T cell activation and cytokine production (Supplementary Fig. 14h). Third, we computed SHAP values, an established interpretability practice for models not employing language embeddings39 and deep neural network learning models37. To ensure robustness, we used three different random seeds (Fig. 7g & Supplementary Fig. 15a) and three data split proportions (Fig. 7h & Supplementary Fig. 15b). In model 3, genes such as tyrosine kinase binding protein (TYROBP), C-C motif chemokine ligand 5 (CCL5), and TCR variable genes (e.g. TRBV13) consistently in the top 10 within CD8.Temra, CD8.GZMK+ Tem and CD4.Temra (Fig. 7g, h & Supplementary Fig. 15). TYROBP is expressed by lymphocytes and mediates signaling for multiple cell-surface receptors40. In addition, TYROBP is highly expressed in certain CD8 T cell subclusters, which are reported to be activated in autoimmune disease41,42. TYROBP has also been linked to M2 tumor-associated macrophage polarization in pancreatic cancer43. CCL5, which induces leukocyte activation and chemotaxis44, has been reported to be upregulated in myocarditis-associated CD8 T cells45, colitis-associated CD4 T cells46 and severe irAE-associated CD4 T cells30. In model 2, TYROBP also rank top in some cell clusters, such as CD8.Temra and CD8.GZMK+ Tem (Fig. 7g, h and Supplementary Fig. 15). Taken together, these results suggest that CD8.Temra likely play a meaningful role in underlying irAE mechanisms.

Therefore, we employed FACS to quantify CD8.Temra within the CD8+ T cell population of PBMCs from two independent ICI patient cohorts (Supplementary Table 5). The first cohort consisted of patients from our scRNA-seq dataset, from whom we thawed freshly frozen PBMC samples. The second cohort included patients with various cancer types undergoing ICI therapy, from whom fresh blood samples were collected. CD8.Temra cells were defined as CD45RA+ CCR7− within the CD45+ CD3+ CD8+ T cell gate (Fig. 8a). Our analysis revealed significantly higher proportions of CD8.Temra in the irAE patients across both cohorts (Fig. 8b–e). These results confirm the association between irAEs and increased CD8.Temra cells are consistent with our cell-type and gene-level importance analysis. To facilitate use of irAE models, we have developed a user-friendly data portal, irAE-atlas (https://iraeatlas.ipm-gba.org.cn). This website not only allows users to calculate the predicted irAE scores for their own samples but also serves as one of the most comprehensive single-cell transcriptome data resources for irAEs.

Methods
Patient populations
In the radio group, patients received weekly intravenous cisplatin for 5 weeks during EBRT, with or without subsequent chemotherapy, depending on physician recommendations. According to the clinical trial protocol (NCT06378840, an observational standard of care treatment study), 17 patients in the radio+ICI group received intravenous nab-paclitaxel and cisplatin weekly for 5 weeks during EBRT, followed by sintilimab 200 mg every 3 weeks for a total of 8 cycles. Participants received anti-PD-1 therapy at no cost as part of their clinical treatment. Three additional patients who were not included in the clinical trial but also received ICI (zimberelimab or penpulimab) after radical chemoradiotherapy were included in the radio+ICI group. Radical radiotherapy consisted of 45 Gy external beam radiotherapy (EBRT) delivered in 5 fractions per week, followed by image-guided brachytherapy which consisted of 30–35 Gy in 5–6 fractions.

We monitored their irAE development of all patients throughout their treatment. We determined the immune relatedness of adverse events by clinical symptoms, laboratory tests and image examinations as well as alleviation of manifestations after discontinuation of ICI. Gastrointestinal irAEs developed following ICI therapy and more than one month after chemoradiotherapy. Laboratory tests were employed to rule out potential infections. Discontinuation of ICI and administration of corticosteroid effectively alleviated the symptoms. For the skin, endocrine, musculoskeletal and connective tissue and hematologic irAEs, they did not occur during chemoradiotherapy, and most were relieved following corticosteroid treatment (Supplementary Table 6). The grading of irAEs was determined for each patient based on the National Cancer Institute Common Terminology Criteria for Adverse Events (CTCAE), version 5. All the anatomical sites where symptoms of irAEs were evident were recorded (Supplementary Table 1).

Single-cell RNA sequencing
For cervical cancer tissues, fresh tissue samples were dissociated at 37°C with gently shaking (50 rpm) using 0.25% trypsin and 10 μg/mL DNase I in PBS containing 5% fetal bovine serum (FBS; Thermo Fisher, Cat. no. SV30087.02) for approximately 40 min. The resulting cell suspension was filtered through a 40 μm nylon strainer to remove debris and clumps. After removing red blood cells by 1× Red Blood Cell Lysis Solution (Thermo Fisher, Cat. no. 00-4333-57), the dissociated cells were then washed and resuspended in 1× DPBS with 2% FBS. For PBMCs, we diluted 2 mL of blood with an equal volume of DPBS and carefully layered this onto 4 mL of lymphocyte separation medium in a 15 mL tube. After centrifugation at 800 g for 20 min, we transferred the mononuclear cell layer to a new tube and added 5 mL of DPBS with 2% FBS. The cells were then centrifuged, treated with red blood cell lysis solution, washed and resuspended in 1× DPBS with 2% FBS. We assessed the cell viability by using 0.4% trypan blue (Thermo Fisher, Cat. no. 14190144) staining and a Countess® II Automated Cell Counter (Thermo Fisher). The cell suspension was then processed using the 10× Genomics Chromium controller to generate gel beads in emulsion (GEMs) and this produced barcoded cDNA by following the manufacturer’s instruction (CG000206 RevD). Sequencing libraries were quantified using a High Sensitivity DNA Chip (Agilent) on a Bioanalyzer 2100 and the Qubit High Sensitivity DNA Assay (Thermo Fisher Scientific). Libraries were sequenced on a NovaSeq 6000 (Illumina) platform using 2×150 chemistry.

Single-cell 5’ and V(D)J enriched library preparation
We prepared single-cell V(D)J libraries using the 10× single cell immune profiling solution kit. We utilized the Chromium Single-Cell V(D)J Enrichment Kit to prepare full-length TCR VDJ segments and full-length BCR VDJ segments, adhering to the manufacturer’s guidelines. The process began with washing and resuspension of the cells in PBS containing 0.04% BSA to achieve a concentration of 100–800 cells/μL. We then captured 500–7000 cells in droplets, with estimated multiple rates of 0.4–5.4%. After encapsulating cells for reverse transcription and barcoding, we broke the emulsions formed and purified the cDNA using Dynabeads MyOne Silane. The cDNA underwent PCR amplification under specific conditions. For the 5’ gene expression library, we used 2.4–50 ng of amplified cDNA, which was fragmented, end-repaired and size-selected with SPRI select beads before PCR amplification with sample indexing primers. Similarly, for V(D)J-enriched libraries, we amplified V(D)J transcripts from 2 μL of cDNA, followed by fragmentation, end-repair and size selection. We sequenced single-cell RNA libraries on Illumina’s HiSeq 4000 platforms with a minimum goal of 25,000 reads per cell. For V(D)J libraries, we aimed for at least 5000 reads per cell, and these were sequenced on the Illumina HiSeq 4000 system.

scRNA-seq, scTCR-seq and scBCR-seq data processing
We processed the raw scRNA-seq data using CellRanger software (version 7.1.0). The sequencing reads were aligned to the human reference genome GRCh38 and quantified using default parameters. We then removed ambient RNAs and any technical artifacts by using the CellBender49 ‘remove-background’ function with default parameters and identified and removed the doublets via DoubletFinder50 v2.0.3. To ensure data quality, we applied the following filtering criteria: cells with fewer than 200 expressed genes and those with more than 10% of reads from mitochondrial genes were excluded and removed. After filtering and normalization, our final dataset had an average of 1844 genes and 5586 unique transcripts per cell.

Dimensionality reduction, clustering and annotation
We merged all the samples and applied R package harmony v1.2.0 to remove batch effects. The Seurat (v5.0.1) R toolkit was used to complete dimensionality reduction and clustering. The clusters on a 2D map were further visualized by using uniform manifold approximation and projection (UMAP). We assigned cell type labels to clusters based on canonical marker genes: immune cells (PTPRC), T cells (CD3D), NK cells (NKG7), B cells (CD79A), plasma cells (MZB1), monocytes (FCN1), macrophages (CD68, CD163), neutrophils (CSF3R), dendritic cells (CD1C), mast cells (TPSAB1), fibroblasts (COL1A1), endothelial cells (VWF) and epithelial cells (KRT6A).

Copy number analysis and malignant cell determination
We utilized the R package inferCNV (version 1.2.1), which detects CNVs and provides insights into genomic alterations at the individual cell level, to conduct CNV analysis at the single-cell level. We also used normal epithelial cells as reference cells18 and excluded the sex chromosomes from the analysis. We further distinguished cancer cells from normal epithelial cells within the epithelial cell population by analyzing their cluster distribution, marker gene expression and inferred copy number variations.

scTCR-seq analysis
We utilized the Cell Ranger VDJ pipeline (version 7.1.0) to analyze scRNA-seq data. This process involved assembling T-cell receptor sequences, identifying CDR3 regions and TCR genes, and annotating T-cell clusters. Cells were selected based on their identification as T cells and the presence of both productive TCR α and β chains. Clonal expansion was defined as identical chains from at least three cells. Clonal diversity, abundance and overlap of the TCR repertoire were quantitatively measured using the clonalDiversity, quantContig and clonalOverlap functions in the scRepertoire package51 (version 1.12.0), respectively.

scBCR-seq analysis
To assemble BCR sequences and identify CDR3 sequences and BCR genes, we used the CellRanger VDJ pipeline (version 7.1.0). We selected cells based on their classification as B cells in the scRNA-seq data and the presence of functional BCR heavy and light chains. Clonal expansion was defined as identical BCR heavy and light chains present in at least three cells. In order to analyze SHM, we reconstructed germline sequences of BCRs using CreateGermlines.py function and calculated number of mutations in each sequence using observedMutations function from Shazam v1.1.052. We assessed BCR clonal abundances using the estimateAbundance functions from the alakazam (version 1.3.0) package52.

Xenium tissue preparation and RNA assay
We performed the Xenium in situ assay using the 10× Genomics Xenium platform according to the manufacturer’s guidelines. We prepared 5 µm sections of formalin-fixed, paraffin-embedded tissue microarrays on Xenium slides. The tissues were de-paraffinized and permeabilized to expose mRNA. We then targeted mRNAs using 380 specific probes and 2 controls: one for non-specific binding and another for genomic DNA to ensure the signals originated from mRNA. We incubated tissue slides overnight at 50°C using a 10 nM probe solution. After washing, we ligated the probes for two hours at 37°C and annealed a rolling circle amplification (RCA) primer. We then amplified the circularized probes to create multiple copies of gene-specific barcodes for each RNA-binding event. This process involved incubation for one hour at 4°C, followed by two hours at 37°C. Following amplification, we washed the slides and chemically quenched background fluorescence to reduce autofluorescence from tissue components. The prepared slides were then placed into imaging cassettes for analysis using a Xenium analyzer.

Xenium data processing
The Xenium analyzer captured multi-layered (Z-stack) images for each cycle across all channels. We processed and combined these images to construct a comprehensive spatial map of transcripts in the tissue section. Using the DAPI (a DNA stain) image as a reference, we integrated image stacks from various fields of view and color channels to produce a single, representative image of the entire tissue section. We then used Squidpy (v1.6.0) and Seurat (v5.0.1) software packages to analyze and visualize the Xenium results, respectively. Cells with at least one transcript expressed were kept for subsequent analysis. We assigned cell type labels to clusters based on the canonical marker genes used in scRNA-seq analysis.

Generation of irAE signatures based on DEGs
We adopted methods from Yan et al., who generated an activate TLS signature from scRNA-seq data53 and Xue et al., who developed TIMELASER tumor subtypes from scRNA-seq data54. Consistent with the process used for developing deep learning models, we first identified irAE-associated cell subclusters in blood samples collected at pre-radio, pre-imm and on-imm time points using two criteria: 1) a significantly increased proportion in irAE patients, with FDR < 0.1 as the cutoff, and 2) a total percentage exceeding 5% in the tissue at the specified time points. We then performed DEG analysis between irAE and non-irAE groups within all irAE-associated cell subclusters, retaining genes with a fold change > 2 and adjusted P value < 0.01. To minimize unexpected noise and expression artefacts from single-cell dissociation, we excluded genes associated with mitochondria, heat-shock protein, ribosome and dissociation54. This process yields three lists of upregulated genes in irAE-enriched immune cells from pre-radio, pre-imm, and on-imm blood samples, respectively, referred as signatures 1, 2, and 3.

Generation of irAE predictive models
Given the majority of irAE-related alterations were observed in peripheral blood samples and the diverse molecular features across different cell types, we employed deep learning strategies to develop predictive models for irAEs based on the scRNA-seq data from pre-radio, pre-imm and on-imm blood samples. The model development process involved four key steps: 1) Cell classification. Immune cell subclusters were categorized as irAE-associated (P1) or non-irAE associated (P0) cells based on their distributions. Subclusters significantly enriched in the irAE group were classified P1, while those with lower proportions were classified P0. Statistical significance was determined using the two-sided Dirichlet-multinomial regression model, with FDR values calculated by the Benjamini–Hochberg method (FDR < 0.1 was used as a significant cutoff). Subclusters comprising less than 5% of the tissue were excluded. 2) Gene ranking and sentence formation. For each cell in these two subgroups, non-zero genes were ranked by raw expression counts, and the top 260 genes were selected to form a ‘sentence’ based on their gene names and their ranks. The choice of 260 genes was made through multiple rounds of testing, as excessively long lists added noise, while shorter lists risked omitting critical information. Each sentence s for cell c was represented as f g1, g2, g3, ..., gn g, where n is gene rank (max=260 in this case) and gn denotes the name of n-th genes. The input data for BERT model was I = f(s1, y1), (s2, y0)..., (si, yi)g, where y1 indicated a cell from P1 and y0 indicated a cell from P0. Using these ‘sentences’ for cells or samples allowed better generalization across datasets generated from different platforms. 3) Prototype learning. The embedding function fφ was used to train BERT model, which constructs a network of gene-to-gene relationships. A basic form of prototypical networks55, a meta-learning approach introduced by Jake Snell, was implemented to fairly learn from various cell subclusters in P1 and P0. For each class (k=0 for P0 and k=1 for P1), we compute an M-dimensional representation xk 2 RM through the embedding function fφ : RD -> RM with learnable parameters φ.

4) IrAE score calculation. Predicted irAE scores were calculated by generating gene sentences sn for cell n and applying the embedding function fφ to obtain embedding vectors en = fφ(sn). Using the Euclidean distance function d, we computed r1 = d(en, x1) and r0 = d(en, x0), representing the distances to irAE (P1) and non-irAE (P0). The relative distance r1 - r0 was used to assess irAE risk, with smaller relative distances indicating higher irAE risks.

Predicted irAE scores (ei) for cell i were normalized to a 0–1 scaling using: ei = 1 - (di - min(d)) / (max(d) - min(d))
where di is the relative distance of cell i and min(d) and max(d) are the minimum and maximum relative distances in the entire dataset, respectively.

The BERT model architecture was implemented using the Transformers (v4.47.1) package by Hugging Face. Prototypical networks was built using PyTorch (v2.5.1) in python (v3.11.11). PyTorch was employed for efficient gene ordering and sentence generation.

For the public scRNA-seq datasets (Supplementary Table 2), we used the cell type annotations when the authors provide these. When the cell type annotation was not provided, we used the Python package CellTypist (v1.6.1) to annotate the cells. We calculated the relative distance for each individual cell and normalized them within each dataset.

Gene and subcluster importance estimation
To evaluate gene contributions to predictions of irAEs, we calculated the SHAP value using the Explainer function of the SHAP package (version 0.48.0)56. After initializing our prototype network with pre-trained weights, we implemented a gene-centric tokenization approach, treating each gene symbol as an indivisible token. SHAP values were then rescaled from 0 to 1. For each gene, we averaged the SHAP values across cells within each subcluster and ranked the genes according to their averaged SHAP values. An ablation analysis was performed by deleting cell subclusters one at a time from each irAE predictive model and comparing four metrics (F1 score, ACC, AUC, and PRAUC) to evaluate the impact of each cell type’s removal. AUC, ACC, PRAUC and F1 score were calculated by R package MLmetrics (version 1.1.3).

Cumulative incidence analysis for irAEs
We acquired the bulk RNA-seq data and the corresponding clinical information for patients with cancer who received ICI therapy from publicly available datasets in the Gene Expression Omnibus (GEO) and Sequence Read Archive (SRA) repositories (accession numbers GSE186143 and SRP353971, respectively). Predicted irAE scores were calculated as described in previous steps for each sample by the specified models. We then used the Kaplan-Meier analysis to assess the time to onset of irAEs in the specified patient groups. P values were determined by log-rank test for cumulative incidence analysis. Patients in SRP353971 were divided into high and low groups based on the optimal cutoff derived from the Youden index of the ROC curve. Due to the limited number of patients in GSE186143, they were divided into high and low groups based on the median values.

Bulk RNA-seq
Blood samples were collected from 14 LACC patients after chemoradiation but before the administration of ICIs for bulk RNA-seq. Total RNA was extracted using TRIzol reagent (Thermo Fisher Scientific, Cat. no. 15596026) and sequenced on an Illumina NovaSeq 6000 platform, with 150 bp paired-end reads, by Shanghai Biotechnology Corporation in accordance with the manufacturer’s guidelines. The clean data generated from the FastQ raw data were subsequently aligned to the human GRCh38 reference genome using Hisat2 software (version 2.0.4). Gene expression levels were quantified using featureCounts and transcript per million (TPM) values were calculated with the calculateTPM function from the R package scuttle.

Multiplex immunofluorescence (mIF) staining
LACC tissue samples were fixed in 4% paraformaldehyde, embedded in paraffin, and sectioned to a thickness of 4 μm. After deparaffinization, these sections were rehydrated using graded alcohols and transferred to the antigen retrieval solution. Following this, the sections were blocked with 5% BSA and then incubated at room temperature for 40 min with the following primary antibodies: CD3 (1:200, Abcam, ab16669), CD4 (1:200, Abcam, ab133616), CD8 (1:200, Abcam, ab237709), CD20 (1:200, Abcam, ab64088), CD68 (1:400, Cell Signaling Technology, 97778), pan-CK (1:200, Abcam, ab234297). After washing three times with PBS-T over 15 min, the sections were exposed to the corresponding secondary antibody for 20 min at room temperature. Subsequent labeling was performed using the Opal 6-Plex Detection Kits (AKOYA Biosciences, NEL821001KT) according to the manufacturer’s protocols. The same procedure was repeated for the additional antibodies. Nuclear staining was achieved using DAPI Solution (Thermo Fisher Scientific, 62248) with a 10-minute dark incubation at room temperature. Images were captured using a 3DHISTECH scanner (Hungary), and intercellular distances were analyzed using HALO software (V3.6.4134.263).

Flow cytometry
Peripheral blood mononuclear cells (PBMCs) were resuspended in FACS buffer (PBS with 2% fetal bovine serum (FBS)) and centrifuged at 300 × g for 5 min at 4°C. After supernatant removal, 100 μL of L/D stain (0.1%, BioLegend [423101]) was added to each tube and incubated for 15 min at room temperature. Cells were then washed with 0.5 mL PBS, centrifuged at 350 × g (5 min, 4°C), and treated with 50 μL of human FcR blocking reagent (0.4%) for 15 min at room temperature to reduce nonspecific binding. Surface antibodies (BioLegend: hCD45RA [304127], hCCR7 [353207], hCD8 [344739], hCD3 [300308], hCD45 [368504]) were diluted 1:100 in 50 μL staining solution, added to cells, and incubated for 30 min at 4°C in the dark. After incubation, cells were washed with 0.5 mL FACS buffer, centrifuged (350 × g, 5 min, 4°C), resuspended in 0.5 mL PBS, and filtered. Cell populations were acquired using a BD FACSAria™ III flow cytometer.

Statistical Analysis
We performed all statistical analyses utilizing R (version 4.3.0). We employed a Dirichlet-multinomial regression model57 to compare cell percentages within each sample. Because the total percentages of all cell subcluster equals 100, an increase in the percentage of one cell subcluster will inevitably cause a decrease in the percentages of other cell subclusters. The Dirichlet-multinomial regression model evaluates differences in cellular distribution between irAE groups (e.g., irAE patients versus non-irAE patients) while considering the proportions of all of the other cell types. The other comparisons between two groups were made with the two-sided Wilcoxon rank-sum test unless otherwise stated, and program enrichment significance among subclusters was determined using the two-sided Kruskal-Wallis test. For cumulative incidence analysis, we employed the Kaplan-Meier method and log-rank to identify the statistical difference. We considered FDR values < 0.1 and/or P values < 0.05 to be statistically significant.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The authors have appropriately addressed my concerns. I retain some concern that the findings will not validate consistently, or are not completely uncoupled from the clinical response, but the authors have shown at least a reasonable semblance of data that they may be useful (the models) in future studies.", "reviewer_2": "I appreciate the authors' comprehensive responses to my initial comments. The manuscript has been revised and improved. However, there are still a few points that require further clarification before I can fully recommend this work for publication:\nMajor Points\n1. Model Interpretability and robustness:\nWhile the authors have made progress in explaining their approach, the model's interpretability remains limited. The transformer-based approach makes it difficult to pinpoint exactly which biological mechanisms are driving the predictions. I am concerned about the generalizability of the models and whether they're capturing dataset-specific rather than biologically meaningful signals.\nI understand that DEG-based irAE signatures are limited, but I would like to see how reliable the transformer model’s predictions ultimately are. Can the authors provide additional evidence, such as ablation or subsampling studies, to identify which features contribute to the superior performance and whether the results are robust.\nThe following statement seems too strong and may rather be a question of sample size and sample complexity, so it should be adapted or removed:\n“These results indicate that DEG-based irAE signatures cannot fully capture the molecular complexity for irAEs, further highlighting the intricate nature of immune responses in irAEs. This complexity also justifies our decision to use deep learning strategies”\n2. Significance Scores for Cell Subclusters:\nThe feature deletion analysis to determine significance scores for cell subclusters is informative, but raises additional questions:\n-- How stable are these significance scores across different model initializations?\n-- Does the relative importance of subclusters change when using different random seeds or data splits?\nIn summary, while the authors have addressed many of my initial concerns, further clarification on model robustness, generalizability, and interpretability would strengthen the manuscript's impact and reliability.", "reviewer_3": "Thank you for the opportunity to review this revised manuscript. I am impressed by the comprehensive single-cell, spatial, and V(D)J profiling study of patients with locally advanced cervical cancer receiving chemoradiation and immune checkpoint inhibition, with a focus on understanding the immune mechanisms underlying immune-related adverse events (irAEs). The scope of this work is impressive, and the longitudinal design with paired tumor and peripheral blood samples represents a valuable resource.\nHowever, I share concerns expressed during the initial review regarding the potential for overinterpretation of the data given the relatively small cohort size and the complexity of the immune dynamics being studied. While the analyses are technically strong, several key issues remain that need to be addressed before the conclusions can be fully supported.\nMajor Concerns\nClarify Attribution of irAEs vs. Radiation-Related Toxicities\nGiven the known overlapping toxicity profiles of radiation and immunotherapy (especially in pelvic malignancies), it is essential to clearly describe how the authors distinguished irAEs from radiation-induced side effects. The manuscript states that events were deemed immune-related based on clinical symptoms and response to ICI discontinuation, but this approach is subjective and lacks detail.\nPlease clarify whether organ-specific attribution, imaging, biopsy, biomarker, or temporal patterns were used to separate radiation-related effects from irAEs. This distinction is foundational to the validity of the downstream analyses.\nSeparate Analyses of irAEs and Treatment Response\nThe manuscript at times conflates immune signatures associated with irAEs and those associated with clinical outcomes such as relapse or progression. While both are important, they are biologically distinct and likely driven by different immune mechanisms.\nI recommend separating these analyses more explicitly and avoiding interpretation that links immune cell subsets or transcriptomic features to both efficacy and toxicity without direct evidence.\nFocus and Validate a Key Immune Signature\nThe dataset is rich, but the analyses attempt to draw conclusions from a very large number of immune cell subclusters and features in a small patient cohort. This raises the risk of false positives and overfitting.\nI suggest focusing the conclusions around one or two key findings (e.g., CD8+ Temra expansion, HPV-specific TCRs) and attempting validation in a second cohort or through an orthogonal method (e.g., flow cytometry, IHC).\nThe deep learning model is intriguing, but its interpretability is limited. Consider including gene-level or subcluster-level feature importance to clarify what drives the prediction.\nStratify irAEs by Organ System or Grade\nGrouping all irAEs together assumes a shared mechanism, which is unlikely. Endocrinopathies, dermatitis, and colitis, for example, have distinct pathophysiologies.\nWhile power is limited, an exploratory analysis stratifying by irAE subtype (or at least separating cutaneous vs. visceral vs. endocrine) could help refine the immune correlates and improve the specificity of predictive signatures.\nAdditional Suggestions\nThe language in some sections overstates the strength of associations (e.g., “underscoring,” “highlighting the role of”).\nPlease consider tempering conclusions and using more cautious phrasing such as “associated with” or “may suggest.”\nThe online tool for irAE prediction is a strength of the manuscript and has translational potential, but further clarity on how clinicians or researchers should interpret or act on these predictions would improve the utility of this resource.\nSummary\nIn summary, this study has great potential, and the authors should be commended for the technical rigor and thoughtful design. However, given the complexity of the biology and the limitations of sample size, I recommend a more focused, hypothesis-driven approach that avoids conflating irAE and treatment response analyses. Greater attention to attribution of irAEs and efforts to stratify or validate core findings would greatly strengthen the manuscript."}, "author_response": "Response: We sincerely appreciate the recognition of the effort invested in data collection of our study, particularly as it serves as a resource paper.\n1. The sample size is very small, and includes only 20 patients treated with ICI. The study also assumes that all irAEs are 'equal' in their mechanism, and we know this is likely not the case from a mechanistic body of literature. The effect sizes are very small, and not corrected for FDR, and in many cases driven by 1 or 2 data points (see figure 2e for example) and are likely therefore not very stable.\nResponse: We appreciate the reviewer’s comments.\nIn this study, we conducted a prospective study (NCT06378840) to recruit patients and closely monitored them throughout the entire treatment process. Given the heterogeneity of irAEs in both time of onset and manifestations1-5, it remains challenge to accurately determine irAEs. Our research group has successfully collaborated with clinicians to properly identify patients with irAEs (Nature Communications 2020; Journal for Immunotherapy of Cancer 2023) and has discussed the challenges associated with sample collection in irAE studies (Nature Reviews Clinical Oncology 2022). The prospective design and stringent protocols of this study enabled us to reliably identify irAEs, which, as you noted, required a significant investment of effort. Consequently, obtaining a large number of patients with fresh samples for such a study is particularly challenging.\nIt is worth noting that recently published high-impact irAE studies employing scRNA-seq strategies have collected tissues, many of which are frozen samples, from tens of patients. Our sample size is comparable to these works (Table R1). Most importantly, our study collected 96 fresh tumor and blood samples for scRNA-seq, which may yield higher cell integrity and RNA quality6. Our sampling strategy provides the advantage of generating data capable of uncovering irAE-associated biological signals that may be missed in studies using frozen samples. Moreover, the prospective nature of our study does not allow for pooled PBMC scRNA-seq via sample multiplexing, as employed by three studies in Table R1, which also requires additional effort and costs. Finally, to the best of our knowledge, our study represents one of the most extensive longitudinal collections of single-cell multi-omics data from patients with or without irAEs, aimed at characterizing molecular features associated with irAEs.\nTable R1. Overview of samples in recently published high-impact irAE studies.\nStudy Journal PubM Sample Pooled ICI patients Non-ICI Cell\nyear ed ID type PBMC (scRNA-seq) patients enrichme\n(Fresh/Fro scRNA-seq (scRNA-seq) nt\nzen)\nGupta Cancer 38744 Frozen Yes 37 22 Yes,\net al. Cell, 246 CD45/CD\n2024 3\nBlum Nature, 37790 Some are Yes 55 (13 with 6 Yes,\net al. 2024 460 fresh, some both PBMC CD45\nare frozen and tissue)\nThoma Nature 38724 Some are Yes 21 8 Yes,\ns et al. Medicin 705 fresh, some CD45\ne, 2024 are frozen\nDimitr Nature 39210 Frozen No 6 0 No\niou et Cancer, 005\nal. 2024\nWe agree that irAEs at different anatomic sites may involve different molecular pathways. However, aberrant activation of T cells has been reported as a contributing factor in the development of various types of irAEs7,8. For example, several recently published studies have highlighted the important role of circulating T cells in colitis (Thomas et.al. Nature Medicine 2024)9, myocarditis (Blum et.al. Nature 2024)10 and any irAEs (Dimitriou et.al. Nature Cancer 2024; Lozano et.al. Nature Medicine 2022)11,12. Given this, we hypothesized that there may be common immune response pathways in immune cells that contribute to any irAEs5. To test this, our study collected tumor and blood samples, not irAE affected tissues, at different treatment stages to comprehensively characterize the dynamic features of various immune cells in the TME and circulating system associated with irAEs.\nTo address the reviewer’s concern that a few of our findings may not be stable, we compared our findings about immune features for irAEs with previously published studies and found them to be consistent. For instance, the role of CD8+ effector T cells has also been observed in disease pathogenesis of ICI induced arthritis13 and myocarditis14,15. Similarly, the role of CD4+ effector T cells has been reported in thyroid-irAEs16 and any irAEs11. A decreased proportion of monocytes has been observed in checkpoint inhibitor-related hepatitis17. To further validate the robustness of our predictive models, which were one of the major findings of this study, we collected an in-house cohort (n=14) with pre-imm blood samples from patients treated with combined radiotherapy and ICI. We also found significantly higher levels of predicted irAE scores in patients with irAEs (Figure R1). We included these results as panels c & d in revised Fig. 7.\nFigure R1. Boxplots comparing predicted irAE scores\nfrom model 3 (developed using scRNA-seq data from\nblood samples collected after chemoradiotherapy and\nduring ICI treatment, left) and model 2 (developed using\nscRNA-seq data from blood samples collected after\nchemoradiotherapy but before ICI treatment, right).\nPredicted irAE scores were calculated using bulk RNA-\nseq data from pre-imm blood samples (collected after\nchemoradiotherapy but before ICI treatment) of irAE\nand non-irAE patients in the in-house cohort. P values\nwere determined using the two-sided Wilcoxon rank-\nsum test, and FDR values were computed using the\nBenjamini–Hochberg method.\nFollowing the reviewer’s suggestion, we have added FDR corrections into our statistical analyses and employed an FDR < 0.1 as the significant cutoff (Revised Fig. 1c-d; Fig. 2f; Fig. 3a, c, d; Fig. 4b, c, d; Fig. 5e, f; Fig. 7. Revised Supplementary Fig. 2b; Supplementary Fig. 3; Supplementary Fig. 5f; Supplementary Figs. 6e, f, h, i; Supplementary Figs. 7c, d; Supplementary Figs. 10f, g, h; Supplementary Figs. 13-14). Importantly, our major conclusions remain unchanged, including the observed molecular dynamics of key immune cells and the predictive potentials of irAE models validated across ten datasets.\nCollectively, our study utilized a prospective design to collect fresh tumor and blood samples, enabling accurate identification of irAEs and ensuring higher cell integrity and RNA quality. The sample size is comparable to other high-impact irAE studies, representing one of the largest longitudinal datasets for irAEs. Findings regarding CD8+ and CD4+ effector T cells and monocyte dynamics are consistent with prior studies. Validation with a new in-house cohort further confirmed the superior performance of the irAE models. We have expanded our discussion about the limited sample size (Lines 458-466), copied below for convenience:\n“Secondly, although we analyzed 96 samples at the single-cell level, the sample size of patients remains relatively limited. However, it may require several years of multicenter efforts to prospectively collect longitudinal samples from a large number of patients with irAEs. But the number of irAEs is increasing exponentially in recent years, and potential severity and fatality of irAEs has been realized by the research community recently4,9. Therefore, we employed deep learning methods to efficiently extract representative irAE-related features for irAE predictive models and validated them using 10 datasets. Further studies with a larger sample size are necessary, particular for predicting different types of irAEs.”\n2. Their is no consideration for clinical response, which is known to in general correlate with iRAE. Thus, the very small difference in the final figure even across datasets in predicting irAE outcome isnt clinically useful or meaningful becuase both groups get many irAEs. Not knowing clinical outcome here doesnt make it clear how this information could actually be used should extensive replication and validation occur. More than likely, the mechanistic work stemming from some of the findings could be a more impactful direction, but that would be a very different manuscript.\nResponse: We appreciate the reviewer’s insightful comment. We agree that incorporating clinical responses could enhance the clinical significance of our study. To address this, we expanded our analysis to clinical response in both our single-cell analysis cohort and validation cohort. For patients in our cohort with locally advanced cervical cancer, a good clinical response was defined as the absence of disease relapse two years after receiving first therapy. In validation cohort, we used the durable clinical benefit label provided in the original studies.\nFigure R2. Association between key immune cells and irAEs or disease relapse. a-b, Heatmaps illustrating the association between the percentage of major cell clusters and irAEs (a) or disease relapse (b). The magenta and cyan squares indicate the higher and lower percentages in irAE patients, respectively. The red and blue dots indicate the higher and lower percentages in non-relapsed patients, respectively. The size of the squares and dots represents statistical significance, with large, medium and small dots indicating FDR of <0.01, 0.01-0.1 and ≥0.1, respectively. X, indicates the data was not applicable. c-d, Proportional analysis of T cell subclusters comparing patients with and without irAEs (c), as well as with and without disease relapse (d) at different time points. e-f, Proportional analysis of B cell subclusters comparing patients with and without irAEs (e), as well as with and without disease relapse (f) at different time points. g-h, Proportional analysis of myeloid cell subclusters comparing patients with and without irAEs (g), as well as with and without disease relapse (h), at different time points. The box plots indicate the median±1 quartile, with whiskers denote 1.5× interquartile range. Each dot represents an individual patient. P values were calculated using the Dirichlet-multinomial regression model and FDR values were calculated by the Benjamin–Hochberg method.\nAs shown in Figure R2, comparisons between irAE and non-irAE patients revealed elevated proportions of B cells, T cells and macrophages in blood samples collected from irAE patients prior to ICI administration. After ICI treatment, T cells increased, while macrophages decreased, in irAE patients (Fig. R2a). No significant differences in immune cell proportions were observed in tumor samples at pre-radio and on-radio time points between the irAE and non-irAE patients (Fig. R2a). Interestingly, pre-radio tumor samples from patients who remained relapse-free two years after receiving first therapy showed increased proportions of several immune cells, including B cells, T cells and macrophages. In contrast, blood samples from non-relapsed patients collected prior to ICI treatment displayed elevated levels of B cells and T cells, with macrophages decreased before and after ICI treatment (Fig. R2b). These results suggest that immune responses within the tumor microenvironment may be more strongly associated with treatment efficacies, while systemic immune responses appear to be more closely linked to irAEs. Consistent with the reported role of memory CD4+ T cells and cytotoxicity CD8+ T cells in irAEs11,18, the proportions of CD4.Temra and CD8.Temra, with cytotoxicity and memory phenotypes, were elevated in patients experiencing irAEs, indicating the involvement of both CD4+ and CD8+ T cells in irAE development. However, only CD4+ T cell subclusters, CD4.Terma, CD4.CAPG+ Tm and CD4.TNFRSF9+ Treg, with cytotoxicity and exhaustion signature expression were increased in non-relapse patient (Fig. R2c-d). Unlike T cells, B cell subclusters exhibited fewer significant frequency changes across patient groups. Proportions of Bn-TCL1A significantly increased following ICI treatment, especially in irAE patients. In contrast, memory B cells were more abundant in pre-radio tumor samples from patient without disease relapse (Fig. R2e-f). Among monocyte subclusters, Mono.ISG15 and Mono.FOLR3, which expressed the highest inflammatory response signature, were enriched in pre-imm and on-imm blood samples from irAE patients. Two macrophage subclusters with higher phagocytosis-related signature expression, Macro.TPSAB1 and Macro.TREM2, were enriched in baseline tumor samples from non-relapsed patients (Fig. R2g-h). We also performed analysis about clonality of T cell and B cells, isotype of B cells and function of myeloid cells between irAE and disease relapse subgroups (presented in revised Fig. 4c-g, 5f, revised supplementary Fig. 6e, 10h).\nIn summary, we observed a marked increase of CD4.Temra and CD8.Temera cells in irAE patients following chemoradiation and ICI therapies, emphasizing the central role of T cells in the development of irAEs. Interestingly, the effects of chemoradiotherapy and ICIs on B cell and myeloid cells differed significantly in irAE patients. Chemoradiotherapy led to an increased percentage of B cells with IgA/IgG isotypes, reduced BCR clonality and elevated expression of inflammation-related pathways of myeloid cells in irAE patients. In contrast, addition of ICI treatment led to a higher proportion of Bn-TCL1A cells and an increase expression of phagocytosis-related gene signatures in myeloid cells. However, in non-relapse patients, baseline tumor samples revealed an enrichment of CD4+ T cell with cytotoxicity and exhaustion phenotypes, B cell with memory phenotype and macrophages with phagocytic phenotype.\nIn the validation cohort, we found that GSE186143 provided clinical information about clinical benefits and irAEs. Therefore, we examined the relationship between irAE predictive models and the efficacy of ICI treatment. Interestingly, the predicted irAE scores from three models showed no significant differences between patients who experienced durable clinical benefits and those who did not (Figure R3), suggesting unique predictive potential of irAE-related features by these irAE predictive models.\nFigure R3. Comparison of the predicted irAE scores in patients (GSE186143) with different treatment efficacies. Upper panels: Violin plots comparing the predicted irAE scores in blood samples collected before ICI treatment from patients who either demonstrated durable clinical benefits or those with no benefits. Predicted irAE scores were determined by the irAE models 1-3 (developed using scRNA-seq data of pre-radio, pre-imm blood and on-imm blood samples, respectively). P values were determined by the two-sided Wilcoxon rank-sum test. Lower panels: durable clinical benefit-free survival analysis using the Kaplan–Meier method. The median values were used to classify patients into high and low groups.\nOverall, single-cell multi-omics data-based analyses suggest that chemoradiotherapy and ICIs reshape key immune cells within the systemic immune system, potentially leading to cell killing, autoantibody production and inflammatory response that facilitate irAE development. For treatment efficacy, cell-killing-related functions of key immune cells in baseline TME could serve as indicators of improved patient prognosis. Predicted irAE scores, derived from models developed using blood samples, showed no significantly differences between patients with varying clinical benefits. This further suggests that the immune mechanisms driving irAEs may be mechanistically distinct from those responsible for anti-tumor efficacy, aligning with our observations in lung squamous cell carcinoma patients receiving ICIs19. Figure R2 were included as panels in revised Fig. 1d-e, Fig. 2f, Fig. 4b, Fig. 5e and revised Supplementary Fig. 5f. Figure R3 were included as revised Supplementary Fig. 15. These descriptions were incorporated in Result and Discussion section in revised manuscript (pages 6-12, 15).\n3. The validation is nice, but they are entirely different collection points and methods and its not clear how the authors went from a scRNAseq predictive model to a bulk RNAseq model, which has very different implications (i.e. genes present in multiple immune compartments and the composition of those compartments etc). This makes me worried that their signature is likely an inflammatory signature related to both immune response and priming for the tumor as well as irAEs, which returns to the utility question in point #2 above.\nResponse: We appreciate this thoughtful comment. Given that the majority of irAE-related alterations were observed in peripheral blood samples and that different cell types exhibited divergent molecular features, we employed deep learning strategies to develop predictive models for irAEs based on the scRNA-seq data from blood samples collected at three time points: pre-radio (before chemoradiotherapy and ICI, referred to as model 1), pre-imm (after chemoradiotherapy but before ICI, referred to as model 2), and on-imm (after chemoradiotherapy and during ICI treatment, referred to as model 3). To ensure these models accurately captured representative features of irAEs at various treatment stages, we validated them using samples collected at time points comparable to those used for corresponding model development.\nSpecifically, we used blood samples collected prior to ICI (pre-ICI blood) to validate model 1 and 2, blood samples collected during ICI treatment (on-ICI blood) to validate model 3 and irAE affected tissue, where almost all irAE-related immune responses could occur, to validate all models. We added illustrations about our validation strategies, and also included two additionally publicly available datasets and one in-house cohort (Figure R4). We incorporated Figure R4 as revised Fig. 7 and added the corresponding descriptions to the revised manuscript (Lines 367-372).\nFigure R4. Validation of the irAE predictive models using independent scRNA-seq and bulk RNA-seq datasets. a, Violin plots displaying predicted irAE scores in blood samples collected before (pre-ICI blood) or during ICI treatment (on-ICI blood) from irAE or non-irAE patients. Scores were calculated using the irAE models developed from scRNA-seq data of blood samples collected at three time points: pre-radio (before chemoradiotherapy and ICI, referred to as model 1), pre-imm (after chemoradiotherapy but before ICI, referred to as model 2), and on-imm (after chemoradiotherapy and during ICI treatment, referred to as model 3). P values were determined using two-sided Wilcoxon rank-sum test, and FDR values were calculated using the Benjamini–Hochberg method. b, Heatmap showing changes in predicted irAE scores (calculated by models 1-3) between irAE-affected tissues and normal controls. Red and blue squares indicate increased and decreased scores in irAE-affected tissue, respectively. P values were determined using the two-sided Wilcoxon rank-sum test, and FDR values were calculated using the Benjamini–Hochberg method. * represents FDR<0.05. c-d, Boxplots comparing predicted irAE scores from the model 3 (c) and model 2 (d) using bulk RNA-seq data of pre-ICI blood samples from irAE and non-irAE patients across publicly available datasets and an in-house cohort. For GSE186143 and in-house cohort, P values were determined using the two-sided Wilcoxon rank-sum test, while for SRP353971, the one-sided Wilcoxon rank-sum test was applied. FDR values were computed using the Benjamini–Hochberg method. e-f, Cumulative incidence analysis of severe irAE (GSE186143) or any irAE (SRP353971) onset. Patients were stratified based on predicted irAE scores from model 3 (e) and model 2 (f). P values were determined by log-rank test for cumulative incidence analysis. HR: hazard ratio.\nThe high costs and significant efforts required for scRNA-seq often limit the sample size of such studies. Thus, to broaden the application of scRNA-seq-based predictive models, we tried to apply these to bulk RNA-seq data. Several studies have successfully generated signatures from scRNA-seq and applied them to bulk RNA-seq datasets20-22. For example, Yan et al. generated an activate TLS signature from scRNA-seq data of lung cancer patients and evaluated its expression in bulk RNA-seq data, finding that higher expression correlated with better clinical responses and overall survival rates (Nature Genetics, 2024)20. Similarly, Xue et al. applied single-cell-based TIMELASER subtypes to published bulk RNA-seq data, achieving consistent results (Nature 2022)21. To confirm that the gene expression profiles of bulk RNA-seq data align well with scRNA-seq data, we performed bulk RNA-seq on a subset of samples from the scRNA-seq cohort. Using Spearman correlation analysis, we compared gene expression values from samples with both scRNA-seq and bulk RNA-seq data in our cohort and GSE186144/ GSE186143. The results showed a significantly positive correlation between averaged scRNA-seq and bulk RNA-seq gene expression values (Figure R5), indicating similar gene rank orders. Specifically, since only ranks of genes were used in the training and predicting steps of irAE model development, applying scRNA-seq-based models to bulk RNA-seq could produce consistent results.\nFigure R5. Spearman correlation analysis of gene expression between scRNA-seq and bulk RNA-seq data from samples in our cohort (left) and GSE186144/GSE186143 (right). Averaged scRNA-seq expression for each sample was calculated using the AverageExpression function from Seurat R package. For bulk RNA-seq data, TPM values were used to represent gene expression. All samples included in this analysis achieved an FDR<0.05.\nTo further address the concerns that irAE predictive models are likely an inflammatory signature related to both immune response and priming for the tumor, we exported the top 100 genes (determined by key word extraction approach for BERT model23) from model 2 and model 3, which showed better predictive performance in the validation. We found that these genes were enriched in anti-tumor immune pathways, including T cell activation and cytokine production (Figure R6). Meanwhile, we also observed gene enrichment in pathways not explicitly associated with anti-tumor responses, such as response to virus and mononuclear cell differentiation (Figure R6). Most importantly, the predicted irAEs scores from model 2 and 3 showed no significant differences between patients who experienced durable clinical benefits and those who did not (Figure R3). This suggests that distinct immune responses underlie irAE development and anti-tumor response. Although the irAE predictive models included genes involved in inflammatory response, these models are not associated with clinical response. We added Figure R6 as revised Supplementary Fig. 14c and descriptions in Lines 396-399 in revised manuscript.\nFigure R6. Gene ontology biological process (GO BP) enrichment analysis for the top 100 genes in model 3 (left) and model 2 (right).\nI think the data are good and may be useful to the field, but the present impact of the work herein as written is very preliminary, wthout a practical output.\nResponse: In response to reviewers’ comments, we have extensively revised our manuscript. For example, we included comparisons of irAEs and clinical responses to characterize irAE-specific molecular features and disease relapse-specific molecular features. We also collected an additional in-house cohort and two publicly available datasets to validate the models to validate the irAE models. Moreover, we conducted analyses comparing predicted irAE scores in patients with different clinical responses to ensure the specificity of irAE predictions. Additionally, we exported the top genes from our predictive models and found that they were enriched in anti-tumor immune pathways and also in pathways not explicitly associated with anti-tumor responses. To be best of our knowledge, the user-friendly data portal irAEatlas is the first irAE-specific single-cell level data resource, providing a comprehensive molecular landscape of blood, tumor and irAE affected tissues of irAE patients. To enhance the accessibility of irAE model, we incorporated the algorithm into the website. This allows researchers to calculate predicted irAE scores for their own samples within minutes, making the platform especially useful for quickly gaining insights into irAE risks.\nReviewer #2 (Remarks to the Author):\nSummary\nThe authors describe a study of immune cell dynamics during a combination of chemoradiotherapy and immune checkpoint inhibitor therapy in patients with locally advanced cervical cancer (LACC). The data acquired to study the dynamics comprise single-cell RNA sequencing and single-cell V(D)J sequencing data, acquired at different time points prior and post chemotherapy, and prior and post immunotherapy, respectively. The authors then develop a deep learning model to assign a score for the likelihood of immune-related adverse events (irAEs) based on scRNA-seq of peripheral blood. We subsequently focus on the scRNA-seq analysis and methodological aspects of the analysis given our expertise. While irAEs are an important risk factor in cancer treatment and better understanding the factors contributing to it is an important and topical issue in cancer research, there are several major issues that should be addressed before the manuscript could be considered for publication. In particular, the methods section needs much more detail so the reader can follow the analysis steps taken by the authors.\nResponse: We thank the reviewer to acknowledge the importance of irAEs for cancer treatment field. Indeed, our major goal in this study is to provide a valuable resource for irAE research community and develop predictive models for irAEs.\nMajor\n1. The authors should try better to interpret the results they find in the analysis of the irAEs that contribute to the differences in immune response. In particular, what are the underlying genes or pathways that contribute to the differences observed?\nResponse: We thank the reviewer for this helpful suggestion. In response to this comment, we extensively revised the interpretation of differences in immune responses in the manuscript (page 6-12). Some revised details are listed below for convenience:\nLines 264-267:\n“Overall, these findings indicate that chemoradiotherapy reduced B cell clonality while increasing the abundance and SHM frequencies of IgA+ and IgG+ B cells in blood and tumor samples of patients with irAEs. This dynamic may play a role in the formation of autoantibody, ultimately leading to irAEs.”\nLines 306-314:\n“In summary, chemoradiotherapy combined with ICI treatment increased the frequencies of monocyte subclusters with heightened inflammatory responses and enhanced the phagocytic capabilities of myeloid cells in patients with irAEs. In non-relapsed patients, macrophage subclusters associated with phagocytosis were enriched, while those associated with inflammatory phenotypes were reduced in baseline tumor samples. Angiogenesis-related signatures were consistently elevated in non-relapsed patients across all time points during treatment. These findings suggest that the inflammatory status of monocytes may contribute to irAE development, while angiogenic potential of myeloid cells may positively influence treatment efficacy.”\nLines 423-428:\n“These findings suggest that chemoradiotherapy and ICIs reshape key immune cells in the systemic immune system, potentially leading to cell killing, autoantibody production and inflammatory response that contribute to the development of irAEs. For treatment efficacy, cell-killing-related functions of key immune cells in baseline TME could serve as indicators of improved patient prognosis.”\n2. Related, but different, the authors should establish whether a irAR score could be generated using simpler statistical measures or more interpretable methods such as differential gene expression rather than a hard-to-interpret transformer model.\nResponse: We strongly agree with the reviewer’s point and have therefore attempted to develop irAE predictive models based on differential gene expression (DEG) analysis. We adopted methods from Yan et al., who generated an activate TLS signature from scRNA-seq data 20 and Xue et al., who developed TIMELASER tumor subtypes from scRNA-seq data21. In line with the process used for developing deep learning models, we first identified irAE-associated subclusters in blood samples collected at pre-radio (before chemoradiotherapy and ICI therapy), pre-imm (after chemoradiation but before ICI therapy) and on-imm (after chemoradiotherapy and during ICI treatment) time points using two criteria: 1) a significantly increased proportion in irAE patients, with FDR<0.1 as the cutoff, and 2) a percentage exceeding 5% in the tissue at the specified time points. We then performed DEG analysis between irAE and non-irAE groups within all irAE-associated subclusters, retaining genes with a fold change > 2 and adjusted P value < 0.01. To minimize unexpected noise and expression artefacts from single-cell dissociation, we excluded genes associated with mitochondria, heat-shock protein, ribosome and dissociation21. This process yields three lists of upregulated genes in irAE-enriched immune cells from pre-radio, pre-imm and on-imm blood samples, respectively, referred to as signature 1, 2 and 3. We calculated the average expression of all genes in each list to generate irAE scores of signatures 1-3. When validating the performance of these irAE scores, we found that the scores derived from any of the three signatures did not consistently increase in irAE patients across the validation cohort, regardless of whether scRNA-seq or bulk RNA-seq datasets were used (Figure R7). Cumulative incidence analysis suggested that low irAE scores indicated higher irAE risk or failed to achieve statistical significance (Figure R7). These results indicate that DEG-based irAE signatures cannot fully capture the molecular complexity for irAEs, further highlighting the intricate nature of immune responses in irAEs. This complexity also justifies our decision to use deep learning strategies. We have included Figure R7 as revised Supplementary Fig. 13 and added corresponding text to revised manuscript (Lines 348-357).\nFigure R7. Validation of the DEG-based irAE signatures using independent scRNA-seq and bulk RNA-seq datasets. a-b, Heatmap showing changes in irAE scores (calculated by signatures 1-3) of pre-ICI blood and on-ICI blood (a) between irAE and non-irAE patients, irAE-affected tissues and normal controls (b). Scores were calculated using the irAE signatures developed from scRNA-seq data of blood samples collected at three time points: before chemoradiotherapy and ICI therapy (refereed as signature 1), after chemoradiation but before ICI therapy (refereed as signature 2), and after chemoradiotherapy and during ICI treatment (refereed as signature 3). P values were determined using two-sided Wilcoxon rank-sum test, and FDR values were calculated using the Benjamini–Hochberg method. Red and blue squares indicate increased and decreased scores in irAE patients or irAE-affected tissue, respectively. *, FDR<0.05. c-e, Boxplots comparing predicted irAE scores from the signature 3 (c), signature 2 (d) and signature 1 (e) using bulk RNA-seq data of pre-ICI blood samples from irAE and non-irAE patients across publicly available datasets and an in-house cohort. P values were determined using two-sided Wilcoxon rank-sum test, and FDR values were calculated using the Benjamini–Hochberg method. f-g, Cumulative incidence analysis of any irAE (SRP353971, f) or severe irAE (GSE186143, g) onset time. Patients were stratified based on irAE scores from signatures 1-3. P values were determined by Kaplan-Meier analysis. HR: hazard ratio.\n3. Much more details about the data analysis, in particular the way the irAE score was computed, are needed, e.g. how was the distance computed, what form did the prototype networks take? This information is required to properly assess this part of the analysis.\nResponse: We appreciate the careful consideration and suggestion from the reviewer. We have revised the descriptions in Methods section (Lines 637-677) and add an illustration to describe the process of our method (Figure R8). Revised text is listed below for convenience:\n“Given the majority of irAE-related alterations were observed in peripheral blood samples and the diverse molecular features across different cell types, we employed deep learning strategies to develop predictive models for irAEs based on the scRNA-seq data from pre-radio, pre-imm and on-imm blood samples. The model development process involved four key steps: 1) Cell classification. Immune cell subclusters were categorized as irAE-associated (𝑃 ) or non-irAE associated (𝑃 ) cells based on their distributions. Subclusters that were significantly enriched in the irAE group were classified as 𝑃 , while those with lower proportions were classified as 𝑃 . Statistical significance was determined using the Dirichlet-multinomial regression model, with FDR values calculated by the Benjamini–Hochberg method (FDR<0.1 was used as significant cutoff). Subclusters comprising less than 5% of the tissue were excluded. 2) Gene ranking and sentence formation. For each cell in these two subgroups, non-zero genes were ranked by raw expression counts, and the top 260 genes were selected to form a ‘sentence’ based on their gene names and ranks. The choice of 260 genes was through multiple rounds of testing, as excessively long lists added noise, while shorter lists risked omitting critical information. Each sentence 𝑠 for cell 𝑐 was represented as {𝑔 ,𝑔 ,𝑔 …,𝑔 }, where 𝑛 is gene rank (max=260 in this case) and 𝑔 denotes the name of 𝑛 -th genes. The input data for BERT model was 𝐼 = {(𝑠 ,𝑦 ),(𝑠 ,𝑦 )…,(𝑠 ,𝑦 )}, where 𝑦 indicated a cell from 𝑃 and 𝑦 indicated a cell from 𝑃 . Using these ‘sentences’ for cells or samples allowed better generalization across datasets generated from different platforms. 3) Prototype learning. The embedding function 𝑓 was used to train BERT model, which constructs a network of gene-to-gene relationships. A basic form of prototypical networks24, a meta-learning approach introduced by Jake Snell, was implemented to fairly learn from various cell subclusters in 𝑃 and 𝑃 . For each class (𝑘 = 0 for 𝑃 and 𝑘 = 1 for 𝑃 ), we compute an 𝑀-dimensional representation 𝒙 ∈ ℝ) where 𝑘 ∈ {0,1} through the embedding function 𝑓 :ℝ* → ℝ) with learnable parameters 𝜑.\n1\n𝒙 = = 𝑓 (𝒔 )\n( |𝑆 | ' 𝒊\n(\n(-!,/!)∈2\"\nwhere 𝑺 is the set of sentences in class 𝑘.\n(\nThe Euclidean distance function 𝑑 was used to compute distance of a query cell 𝑐 from the prototypes 𝒙 :\n(\n𝑒𝑥𝑝 (−𝑑(𝑓 (𝒄),𝒙 ))\n' (\np (𝑦 = 𝑘|𝒄) =\n' ∑ 𝑒𝑥𝑝 (−𝑑(𝑓 (𝒄),𝒙 ))\n(# ' (#\n4) IrAE score calculation. Predicted irAE scores were calculated by generating gene sentences 𝑠 for cell 𝑛 and applying the embedding function 𝑓 to obtain embedding vectors 𝑒 = 𝑓 (𝑠 ). Using the Euclidean distance function 𝑑, we computed 𝑟1 = 𝑑(𝑒 ,𝑥 ) and 𝑟0 = 𝑑(𝑒 ,𝑥 ), representing the distances to irAE (𝑃 ) and non-irAE (𝑃 ). The relative distance 𝑟1− 𝑟0 was used to assess irAE risk, with smaller relative distances indicating higher irAE risks. Predicted irAE scores (𝑒 ) for cell 𝑖 were normalized to a 0-1 scaling using:\n#\n𝑑 −𝑚𝑖𝑛 (𝑑)\n#\n𝑒 = 1−\n# 𝑚𝑎𝑥(𝑑)−𝑚𝑖𝑛 (𝑑)\nwhere 𝑑 is the relative distance of cell 𝑖 and 𝑚𝑖𝑛 (𝑑) and 𝑚𝑎𝑥 (𝑑) are the minimum and maximum relative distances in the entire dataset, respectively.\nThe BERT model architecture was implemented using Transformers (v4.47.1) package by Hugging Face. Prototypical networks was built using PyTorch (v2.5.1) in python (v3.11.11). PyTorch was employed for efficient gene ordering and sentence generation.”\nFigure R8. Overview of irAE predictive model development.\nCan the authors find any relationship between the intensity/number of irAEs and the score computed? I understand there are only very few patients, but it would be interesting to at least show these data. I assume that there really is a continuum of irAEs in reality?\nResponse: We appreciate the reviewer’s suggestion regarding the intensity and number of irAEs. To address this, we analyzed three validation datasets containing information about number of irAEs and two datasets with information about irAE grades. While we obtained statistically significant results, the dynamics between different irAE number groups are divergent across datasets. In GSE206301, the highest irAE scores were observed in the group with four irAEs, while the lowest irAE scores were in the non-irAE groups, using models 1-3. However, in the other two datasets, predicted irAE scores did not show a consistent increase in the groups with more irAEs (Figure R9a). For irAE grades, predicted irAE scores calculated using model 2 and 3 were highest in patients with grade 4 irAEs and lowest in patients with grade 1 irAEs or without irAEs (Figure R9b). Overall, these results suggest that irAE scores derived from model 2 and 3 may partially reflect the severity of irAEs rather than numbers of irAEs.\nFigure R9. Comparison of predicted irAE scores using models 1-3 between irAE number groups (a) and grade groups (b). P values were determined by Kruskal-Wallis test. Scores were calculated using the irAE models developed from scRNA-seq data of blood samples collected at three time points: pre-radio (before chemoradiotherapy and ICI, referred to as model 1), pre-imm (after chemoradiotherapy but before ICI, referred to as model 2), and on-imm (after chemoradiotherapy and during ICI treatment, referred to as model 3).\n5. The categorization of immune cells into irAE or non-irAE cells needs more detail and a justification of why this assignment of subclusters to the irAE group is valid. Can the authors provide any rationale, e.g. are any of the genes in the subclusters associated with immune responses? How was the significance of clusters determined?\nResponse: Our results derived from single-cell analysis indicate that chemoradiotherapy and ICIs reshape the composition and activity of key immune cells in the tumor microenvironment and systemic immune system, potentially leading to cell killing, autoantibody production and inflammatory response that contribute to the development of irAEs. For example, following chemoradiation and ICI treatment, patients with irAEs exhibited increased frequencies of monocyte subclusters with heightened inflammatory responses and enhanced phagocytic capabilities of myeloid cells. Additionally, the proportions of CD4.Temra and CD8.Temra, which expressed the highest cytotoxicity signatures, were elevated in irAE patients, indicating the involvement of both CD4+ and CD8+ T cells in irAE development. Therefore, to efficiently extract representative irAE-related molecular features, we used immune cells with significantly altered frequencies to build predictive models. We have revised the descriptions about how we categorize immune cells into irAE or non-irAE cells models in Methods section (Lines 640-646). Revised text is listed below for convenience:\n“1) Cell classification. Immune cell subclusters were categorized as irAE-associated (𝑃 ) or non-irAE associated (𝑃 ) cells based on their distributions. Subclusters significantly enriched in the irAE group were classified 𝑃 , while those with lower proportions were classified as 𝑃 . Statistical significance was determined using the Dirichlet-multinomial regression model, with FDR values calculated by the Benjamini–Hochberg method (FDR<0.1 was used as significant cutoff). Subclusters comprising less than 5% of the tissue were excluded.”\nTo leverage this comprehensive dataset, we employed deep learning methods to extract key information from all irAE cells, without pre-determining the significance of individual cell subclusters during model development. However, as understanding the contributions of specific cell subclusters to the predictive models could improve our understanding of the biological relevance behind irAE scores, we developed a method to determine significance of each subcluster. Specifically, we performed a feature deletion analysis by removing one cell subcluster from input data and comparing the model’s performance to that of the full dataset. The decrease in recall values was used as a measure of the significance of the excluded cell subcluster (Figure R10). Consistent with our previous findings, Mono.FOLR3, which expressed the highest inflammatory response signature, was enriched in the on-imm blood samples of irAE patients. It also showed the highest significance score in model 3, which was developed using scRNA-seq data of on-imm blood samples. Although frequencies of CD4.Temra and CD8.Temra cell subclusters were elevated in irAE patients, CD8.Temra, which showed higher proportions than CD4.Temra in on-imm blood samples, demonstrated higher significance scores in model 3. In summary, these significance scores align closely with our observations regarding cell frequencies and gene expression patterns, further validating the robustness of our model developing method.\nFigure R10. Bar plots showing the significance scores of cell subclusters used in the development of models 1-3. Each bar represents the difference in recall value when a specific cell subcluster is removed from the input data compared to the model utilizing full input data. Higher significance scores indicate a greater impact of the corresponding cell subcluster on model performance. The irAE models developed from scRNA-seq data of blood samples collected at three time points: pre-radio (before chemoradiotherapy and ICI, referred to as model 1), pre-imm (after chemoradiotherapy but before ICI, referred to as model 2), and on-imm (after chemoradiotherapy and during ICI treatment, referred to as model 3).\nMinor\n-- Grammar: 187 By using of TCR sequences (clonotypes) as molecular barcodes,\n-- More details about the number of patients in the different conditions\nResponse: We thank the reviewer for this comment and have substantially edited grammar problems in our manuscript. In addition, we have added the number of patients in different conditions in revised Fig. 7.\nReviewer #3 (Remarks to the Author):\nThe stated objective of the study by Zhang and colleagues is to investigate immune related AEs and associated mechanisms using a longitudinal single-cell multi-omics analysis of the tumor microenvironment in patients with locally advanced cervical cancers (LACC) enrolled on a trial (NCT06378840) of patients undergoing sequential chemoradiation plus immune checkpoint inhibitors (ICIs).\nOverall this is a novel study that addresses an important knowledge gap regarding the nature and predictors of irAEs and provides a web-based platform for further development of insights in new datasets. However, in its current form the study has serious limitations that detract from its potential value. Specific issues are outlined below.\nResponse: We appreciate the reviewer for recognizing the novelty of our study.\nCritiques:\n1. The conclusion that “These models have the potential to optimize patient selection and improve treatment outcomes for individuals receiving a combination of chemoradiation and ICI therapies,” (line 90-92 and in the discussion) is overstated for this small study based on 30 patients (only 20 of which received IO) in a single cancer type and supported by their data.\nResponse: We thank reviewer for this thoughtful critique. To be noticed, we have validated our predictive models across 10 datasets, including 9 publicly available datasets and one in-house cohort, which demonstrated that our models have unique power in predict irAEs. We have revised the descriptions in the introduction and discussion section. Some revised texts are listed below for convenience:\nLines 90-92:\n“These models may help identify patients with higher risks of developing irAEs and could potentially assist in optimizing patient selection for those undergoing combined chemoradiation and ICI therapies.”\nLines 428-430:\n“Overall, this study provides valuable insights into the complex immune responses elicited by combination treatment and offers possible directions for optimizing clinical benefit for cancer patients.”\nLines 481-482:\n“We have also developed an online tool that may potentially help clinicians in assessing irAE risks for ICI patients.”\n2. Data presentation needs to be more specific in terms of discussion of TME vs. peripheral blood dynamic changes. Currently the authors go back and forth and figures do not clearly address whether the data presented in from peripheral blood or tissue.\nResponse: We appreciate your thoughtful feedback. We have carefully revised all the figures to clearly indicate the data is derived from peripheral blood or tissue. Please see our revised Figs. 1-5.\n3. All the tissue biopsies appear to be from the cervix. As authors note, during treatment there is a significant decline in number of tumor cells and hence the later biopsy timepoints are largely comprised of stromal and immune cells. What is the biological rational for biopsies of the cervix providing information relevant to development of irAEs across a variety of other organs?\nResponse: Thanks for raising this question. In this study, we collected both fresh tumor and peripheral blood samples with locally advanced cervical cancer patients at four time points across the entire treatment period. Tumor samples were collected at pre-radio (before chemoradiotherapy and ICI) and on-radio (during chemoradiotherapy) time points, while peripheral blood samples were collected from pre-radio, pre-imm (after chemoradiotherapy but before ICI) and on-imm (after chemoradiotherapy and during ICI treatment) time points. This design enabled us to longitudinally compare the divergent molecular features of immune cells within the tumor microenvironment and systemic immune systems, particularly in the context of irAEs. Several studies have reported that, in cancer patients treated by ICIs, myocarditis and pneumonitis may result from cross-reactive T cells recognizing shared antigens between tumor and normal tissues14,25. Additionally, inflammatory arthritis induced by ICIs has been linked to CX3CR1hi CD8+ T cells trafficking from blood to joints13. Building on these findings, we hypothesized that radiotherapy exposure reshapes the TME, which in turn induces the alterations in the systemic immune system and may contribute to irAEs in other organs. Indeed, our findings indicate that chemoradiotherapy and ICIs reshaped the composition and activity of key immune cells in both the tumor microenvironment and systemic immune system. For example, chemoradiotherapy reduced B cell clonality while increasing the abundance and SHM frequencies of IgA+ and IgG+ B cells in blood and tumor samples of patients with irAEs. Following chemoradiation and ICI treatment, patients with irAEs exhibited increased frequencies of monocyte subclusters with heightened inflammatory responses and enhanced phagocytic capabilities of myeloid cells in blood samples. Notably, spatial transcriptomic analysis of tumor tissues using the 10× Genomics Xenium platform suggested that immune cells associated with irAEs, including B cells, CD4+ T cells, and CD8+ T cells, were aggregated in the on-radio tumors of irAE patients. Spatial distribution analysis of immune cells from mIF images revealed increased proximities between immune cells in irAE patients compared to non-irAE patients (Figure R11). In conclusion, our study identified molecular dynamics underlying irAEs by dissecting longitudinal tumor and peripheral blood samples.\n4. The authors should discuss potential biases that may be introduced due to the different number of tissue samples (varying for 1 to 5) included from different patients\nResponse: Thanks for this useful comment. We have discussed the different number of tissue samples in our discussion section (Lines 466-470). Revised text is listed below for convenience:\n“Due to limited follow-up time or unsuitable patient conditions, some patients lacked samples of one or two time points. As a result, certain comparisons, such as those between the on-imm blood sample of relapsed and non-relapsed patients, were based on a small sample size (e.g., only two patients in the relapsed group). Therefore, the conclusions drawn from such comparisons should be validated in future studies.”\n5. The model for prediction of irAEs is problematic (Figure 7). While they may be able to claim differences be different groups (e.g. irAE present/absent) the distributions are largely overlapping, especially for the most relevant timepoints which are pre-radio and pre-Imm (since these are the scenarios in which knowledge about risk of ir-AEs could produce actionable results). It is difficult to view on treatment discriminators and information gained from irAE affected tissues as useful for prediction, given these biospecimens originate from patients who are already on IO and have developed irAE, respectively.\nResponse: We completely agree with the reviewer that the time points before ICI treatment are critical for irAE prediction, as these stages provide valuable opportunities to implement interventions or alternative treatment designs for patient benefits. This rationale also motivated us to apply our models to three datasets of pre-ICI blood with bulk RNA-seq data. In this study, we developed deep learning models to predict irAEs using scRNA-seq data from blood samples collected at three key time points: pre-radio (before chemoradiotherapy and ICI, referred to as model 1), pre-imm (after chemoradiotherapy but before ICI, referred to as model 2), and on-imm (after chemoradiotherapy and during ICI treatment, referred to as model 3). The primary objective was to identify the model that best captured the molecular features of irAEs and owned the highest predictive potential. To evaluate the models, we first validated them using scRNA-seq data collected at time points corresponding to those used for model development. Pre-ICI blood samples were used to validate model 1 and 2, while on-ICI blood samples were used to validate model 3. Additionally, irAE-affected tissue, where most irAE-related immune responses may occur, was used to validate all three models (Figure R4a-b). Following this initial validation, we applied model 2 and model 3, which showed better predictive performance in the scRNA-seq validation, to bulk RNA-seq data of pre-ICI blood samples. In all three datasets, both models predicted significantly higher irAE scores in pre-ICI blood samples from irAE patients, with model 3 showing the most pronounced differences (Figure R4c-d). Cumulative incidence analyses further demonstrated better performance of model 3. By assessing the time to irAE onset or severe irAE onset based on follow-up duration after ICI administration, we found that patients with higher irAE scores predicted by model 3 experienced significantly shorter times to irAE onset, with hazard ratios of 1.50 and 2.90, respectively. In contrast, model 2 did not demonstrate the same predictive potential (Figure R4e-f). In summary, these findings demonstrated that higher irAE scores predicted by model 3 in pre-ICI blood samples could potentially identify patients at higher risk of developing irAEs within a shorter time frame.\n6. The assertion that “immune cells associated with irAEs, including B cell, CD4+ T cells and CD8+ T cells, were aggregated in the on-radio tumors of irAE patients, whereas this aggregation was absent in the on-radio tumors obtained from non-irAE patients (Figs. 6a & c)” (lines 290-3) seems anecdotal and not requires more rigorous statistical analysis.\nResponse: We appreciate this insightful comment. To address this, we conducted multiplex immunofluorescence (mIF) staining on tumor tissues collected during chemoradiotherapy from patient who developed irAEs or those who did not. To ensure consistency with data obtained from Xenium platform, we employed markers for total T cells (CD3), CD4+ T cells (CD4), CD8+ T cells (CD8), B cells (CD20), macrophages (CD68) and malignant cells (pan-CK) (Figure R11a). Using HALO software, we analyzed the spatial distribution of immune cells by measuring the intercellular distances between B cells, CD4+ T cells and CD8+ T cells (Figure R11b). The results indicated that an increased proximity between immune cells in irAE patients compared to non-irAE patients. Furthermore, we quantified the number of CD20+ in close proximity to CD4+, CD8+ and CD3+ cells and observed significantly more such cell proximities in tumor samples from irAE patients (Figure R11c), suggesting increased immune cell aggregations and interactions in the tumor microenvironment. We included these results as panels in revised Fig. 6, and corresponding descriptions in revised manuscript (Lines 333-342).\nFigure R11. Spatial analysis of immune cell distributions in tumor microenvironment after chemoradiotherapy. a, Representative images of multiplex immunofluorescence (mIF) staining of tumor samples after chemoradiotherapy from patients with or without irAEs. The markers visualized include CD3 (cyan, T cell marker), CD4 (yellow, CD4+ T cell marker), CD8 (purple, CD8+ T cell marker), CD20 (red, B cell marker), CD68 (blue, macrophage marker), pan-CK (orange, epithelial marker) and DAPI (blue, nuclei). Scale bar = 100 µm. b, Spatial plots illustrating intercellular distances within tumor tissues, analyzed using HALO software (V3.6.4134.263). Each dot represents an individual cell, with lines connecting cells that are proximal to each other (distance <= 50 µm). c, Bar plot showing the number of designated cells with distance <=50 µm for each patient with (n=3) or without irAEs (n=3). Statistical comparisons between the two groups were performed using a two-sided Student's t-test.\n7. Some grammatical errors should be addressed (e.g. lines 87-88).\nResponse: We have thoroughly checked our manuscript about the grammatical errors.\n8. Figure 2a panel seems to be missing\nResponse: We have revised Fig. 2 and include it as Figure R12 for convenience.\n9. Figure 2d is complicated and cumbersome to decipher in it’s current form (at least in my opinion).\nResponse: We have revised Fig. 2 and the legends accordingly based on this suggestion. Revised fig. 2 is copied below as Figure R12:\nFigure R12. Associations of T cells with irAEs. a-d, The uniform manifold approximation and projection (UMAP) plots of 248,557 high-quality T cells (a) and CD8+ T cell subcluster (b), CD4+ T cell subcluster (c) and NK T cell subcluster (d). Each dot represents individual cell and colored by their cell subclusters. e, Stacked bar plots showing the distribution of the T cell subclusters in different tissue type, treatment, irAE status and disease relapse. f, Proportional analysis of T cell subclusters comparing patients with and without irAEs at different time points. The box plots indicate the median±1 quartile, with whiskers denote 1.5× interquartile range. Each dot represents an individual patient. P values were calculated using the Dirichlet-multinomial regression model and FDR values were calculated by the Benjamini–Hochberg method. Tn, naïve T cell; Tm, memory T cell; NK, natural killer; Tem, effector memory T cells; Th, T helper cell; Tex, exhausted T cell; Temra, terminally differentiated effector memory or effector T cells; ICI, immune checkpoint inhibitor; irAEs, immune-related adverse events; Radio, chemoradiotherapy; UK, unknown; NoICI, not treated by ICI.\nReference"}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have appropriately addressed my concerns. I retain some concern that the findings will not validate consistently, or are not completely uncoupled from the clinical response, but the authors have shown at least a reasonable semblance of data that they may be useful (the models) in future studies.", "reviewer_2": "I appreciate the authors' comprehensive responses to my initial comments. The manuscript has been revised and improved. However, there are still a few points that require further clarification before I can fully recommend this work for publication:", "reviewer_3": "The response to reviewers concerns is very detailed and sufficiently addresses the concerns."}, "author_response": ""}]

Decision letter