Spatial Positioning of Immune Hotspots Reflects the Interplay between B and T Cells in Lung Squamous Cell Carcinoma

Intratumoral immune hotspots beyond tertiary lymphoid structures reflect an immunosuppressive microenvironment, different from peritumoral immune hotspots, warranting further study in the context of immunotherapies.


Introduction
Despite a plethora of research dedicated to the study of non-small cell lung cancer (NSCLC) intratumor heterogeneity, deficiencies in our knowledge of the fundamental changes that occur within the tumor microenvironment have impeded progress in this area. Diverse immune responses have been proposed for different NSCLC subtypes, lung adenocarcinoma (LUAD), and lung squamous cell carcinoma (LUSC) based on the analysis of genome and transcriptome (1). Our recent findings resolve the link between spatial variability of lymphocyte infiltration with tumor immune evasion in LUAD (2). However, it is unclear how diverse immune cell subsets cooperate in a compact niche to drive pro-or antitumor responses, or whether they can be effectively measured in histology samples as part of the clinical routine (3).
We have previously demonstrated the importance of studying tumor spatial architecture for defining the clinical relevance of the immune response in estrogen receptor-positive and -negative breast cancer (4,5). Spatial aggregation, rather than the sheer abundance of lymphocytes, is associated with recurrence-free survival in both breast cancer subtypes (4,5). While this highlights the importance of examining not just cell abundance but also spatial patterns that can inform disease prognosis, important questions remain: (i) whether our observation extends to other cancer types beyond breast cancer, such as the highly immunogenic NSCLC; (ii) how these spatial patterns translate into immune cell functions and functional phenotypes; and (iii) how they collectively drive a microenvironment that promotes or suppresses cancer progression. In this article, we used transcriptomic data and deep learning-based digital pathology image analysis to: (i) discover the clinical relevance of spatial organization of lymphocytes in NSCLC; (ii) identify the transcriptomics underpinning of immune spatial heterogeneity; and (iii) uncover the interplay of immune cell subsets underlying distinct spatial organization patterns.

Patients and cohorts
All NSCLC whole-tumor hematoxylin and eosin (H&E) sections from The Cancer Genome Atlas (TCGA) analyzed in this study were downloaded from the TCGA data portal (n ¼ 743 LUAD and 712 LUSC). These sections correspond to the total of 473 patients in TCGA LUAD and 466 in TCGA LUSC. Overall survival data (i.e., times to death or follow-up and vital status) was available for 935 patients with NSCLC, whereas additional clinical data (stage, age, smoking packyears) were only available for 723 patients with NSCLC. The multivariate Cox validation took into account equally distributed discovery (n ¼ 237 LUAD and 231 LUSC) and validation (n ¼ 236 LUAD and 231 LUSC) cohorts.
The 10 external validation samples with multiplex IHC staining were obtained from the Leicester Archival Thoracic Tumor Investigatory Cohort (LATTICe) study, which consists of 1700 University Hospitals of Leicester Trust patients who underwent surgical treatment with curative intent for primary lung tumors. A consultant thoracic malignancy pathologist visually examined 20 H&E slides corresponding to 20 patients with high immune clustering. All 20 cases were analyzed using the immune hotspot pipeline to rank the top 10 cases. These 10 cases were then dearchived, and three serial sections were prepared for T-cell, H&E, and B-cell staining. The validation cohort was ethically approved by an NHS research ethics committee (ref. 14/EM/1159).

Automated H&E image analysis mapped heterogeneous cell populations
We collected four orthogonal data types to test the accuracy of automated H&E image analysis for single-cell classification in wholesection NSCLC images. CRImage (6) identified a total of 206,272,040 cancer cells, 35,733,365 stromal cells, and 49,953,400 lymphocytes in 1,455 H&E sections from frozen samples of 939 TCGA NSCLC patients. The balanced accuracy as an average of specification and sensitivity of the classifier was 82.13% for cancer cells, 89.16% for stromal cells and 87.81% for lymphocytes, comparing against 4,169 cancer, 1,533 stromal cells and 1,999 lymphocytes annotations by a thoracic malignancy pathologist in randomly sampled images (Fig. 1A). Automated cell percentage was more strongly correlated with molecular gene signatures than the independent scoring provided by TCGA pathologists (Supplementary Fig. S1A and S1C). This included a comparison between cancer cell percentages given by automated image analysis and pathologic scores of H&E, and tumor purity measures from the gene expression-based method ESTIMATE (7) and copy number-based ABSOLUTE (8). For stromal cell estimates, we compared H&E-derived scores with stromal scores from gene expression-based MCP-counter (9) and ESTI-MATE (7). There was no pathologic score of lymphocytes available from the TCGA dataset, instead the lymphocyte percentage estimated by automated image analysis was evaluated against immune signatures inferred by TIMER (10). Automatic quantification of lymphocyte percentages significantly correlated with CD4, CD8 T-cell, and B-cell scores (R ¼ 0.31, 0.17, and 0.24, respectively; P < 0.01; Supplementary Fig. S1B)

Bioinformatic analysis
The gene enrichment analysis for the S intra/immune in LUSC was performed using the "TCGAbiolinks" package in R. The RNAsequencing data (all Level 3 raw read counts, available from the TCGA data portal) were used as input to the transcriptome differential gene expression analysis. The outputs of highly enriched genes clustered under the Gene Ontology biological processes for LUSC and LUAD are shown in Fig. 2A and Supplementary Fig. S2A, respectively. The highly enriched immune genes relevant to the immune response/activation of immune response processes in LUSC included DEFB4A, LTF, FCN3, C7, MS4A1, PAX5, PLA2G1B, PLA2G2A, CHGA/CGA, FGB, BPIL1, DEFA5, CR2, PLUNC, C8B, DMBT1, KRT1. The differential gene expression analysis was performed using the 'DESeq2' R package to rank differentially expressed genes with logFC larger than 1 in high versus low groups for the S intra/immune in LUSC and LUAD, with the optimal cut-off point determined on the entire cohort (Supplementary  Table S1; Fig. 2B; Supplementary Fig. S2B). The B-cell signatures were estimated by seven algorithms Danaher (11), TIMER (10), CIBER-SORT Absolute mode (12), MCP-counter (9), quanTIseq (13), EPIC (14), and xCell (15). The patient-level estimates were calculated as the average across slide-level scores of the patient. Comparisons between S intra/immune high and low patients stratified by the median revealed a consistent upregulation of B-cell infiltrates in the high group (Fig. 2C). The tumor mutation burden (TMB) was computed as somatic mutations per Mb using R library "maftools" (16) for 444 TCGA LUSC tumors with mutation annotation format files available. The clonal neoantigen load was derived from a previous work (1).

Survival analysis of TCGA datasets
The investigated scores (S intra/immune , S intra/cancer , S intra/tissue ) were first dichotomized using Kaplan-Meier curves (R "survminer" package) to find an optimal cut-point on the discovery cohort (n ¼ 237 LUAD and 231 LUSC). The optimal cut-off point was computed with maximally selected log-rank statistics. Next, the dichotomized cut-off point found in discovery was used in a multivariate Cox proportional hazards model (R "survival" package) adjusted for stage, age, and smoking pack-years to further substantiate the survival prognosis. The Kaplan-Meier curve and the Cox model with the same cut-off point were then tested in the validation cohort (n ¼ 236 LUAD and 231 LUSC), which revealed that the S intra/immune is a LUSC-specific prognosticator. The survival analysis was repeatedly performed for four patient groups determined by their smoking history (current reformed smokers for less than 15 years, current reformed smokers for more than 15 years, current smokers and lifelong nonsmokers; n ¼ 457 LUAD and 446 LUSC in total) with the cut-off point of S intra/immune optimized for each category. Lymphocyte percentage and TMB were dichotomized using a cut-off point optimized for the discovery cohort. To evaluate the contingency of survival analysis on the split of cohorts, we randomly split the TCGA LUSC cohort into discovery and validation by a 1:1 ratio for 100 times. In each iteration, the numerical parameters (S intra/immune , age, pack years, and TMB) were dichotomized by the threshold optimized for the discovery cohort, followed by a log-rank test conducted in the validation cohort to determine the statistical difference in the OS of patients stratified by the dichotomized parameters. The stage was assessed as a categorical parameter in the validation cohort.

Identification of hotspots
Hotspot analysis was conducted on H&E slides using the method and parameters as previously proposed (5). Getis-Ord spatial analysis was carried out over 50Â50 mm 2 grids using the following equations: where n is the number of grids; C j is the count of cancer cells or lymphocytes in grid j; W i,j is 1 if grid i and j are neighbors, otherwise 0. Grids sharing a vertex or edge with the analyzed grid j were defined as first-order neighbors. Second-order neighbors included grids sharing a vertex or edge with the first order neighbors and so on. Grids within the fourth order neighborhood were considered neighbors to grid j in the equation (5). P value was calculated on the basis of Z score thresholds (17). Regions with P < 0.05 regarding cancer cells and P ≥ 0.05 regarding lymphocytes were categorized as cancer hotspots (CH). Regions with P < 0.05 regarding lymphocytes and P ≥ 0.05 regarding cancer cells were categorized as peritumoral immune hotspot (IH peri ). And intratumoral immune hotspot (IH intra ) was defined as regions with P < 0.05 for both cancer cells and lymphocytes. The definitions of IH peri and IH intra were described by the following equations.
where P c x and P l x represent the P value quantifying the significance level of enrichment of cancer cells and lymphocytes in grid x, respectively.
To quantify the spatial relationship between immune and cancer aggregates for each tumor, we derived three scores based on different normalization terms for the IH intra (5).
1. To quantify the proportion of IHs interfacing with tumor, the intratumoral immune hotspot score S intra/immune was computed as the proportion of IH intra in the immune cell rich area.
2. To measure the amount of intratumoral immune aggregates in relation to the cancer area, S intra/cancer was computed as the ratio of IH intra to the sum of CH and IH intra .
3. To quantify the colocalization of immune and cancer aggregates in the context of the entire tissue, S intra/tissue was computed as the fraction of IH intra in the whole tissue including both hotspot area and non-hotspot area.

Deep learning pipeline to map cell populations on H&E and IHC images of the external validation samples
We performed cell detection and classification on H&E slides and the IHC T-cell panel of the validation cohort using the SCCNN pipeline (2, 18) pretrained on H&E and IHC images, respectively. Because of variations in the staining protocol, we trained an additional SCCNN-IHC classifier to discriminate CD4 þ FOXP3 À and CD4 þ FOXP3 þ cells from this cohort. The SCCNN pipeline consists of two models to perform cell detection and classification sequentially. The cell detection was performed on 31Â31 pixels image patches to predict pixels with high probabilities of being the cell center. The output is a table of predicted cell coordinates without class label. The table is then used to retrieve single cell images of 51Â51 pixels centering at the predicted cell location. The second step involves classifying nine neighboring locations near the center of the cell. The most frequent class predicted among the nine neighborhoods will be determined as the class of the cell (2,18).
Prior to the model implementation, whole slide images of H&E and IHC samples were rescaled to 0.45-mm/pixel and divided into tiles of 2000Â2000 pixels to match the resolution of training images of the SCCNN. A color normalization step was performed for H&E tiles using an example H&E tile from the training set as reference. The H&E tiles were first converted from RGB to Lab color space. Then a Reinhard stain normalization from the Image processing toolbox in Matlab R2019b (RRID:SCR_001622) was used to transform all the H&E tiles based on the reference image. We did not normalize the color for the IHC slides. Instead, the classifier was trained on single cell patches augmented by contrast and brightness adjustment, which encouraged the model to capture stain variability in the dataset.
For the B-cell panel, we trained a distance regularized dense inception network (DRDIN; ref. 19) for cell detection using masks of annotated cell locations, and an SCCNN-IHC classifier (18) to differentiate between three dominant B-cell subsets (CD20 þ CXCR5 À , CD20 þ CXCR5 þ , and CD79b þ ), and the P40 þ neoplastic cells. Similar to SCCNN detector, the DRDIN model was trained to predict pixels of the cell center (19). However, the input image size was of 224Â224 pixels, allowing the integration of features in a larger spatial context and therefore a better capability to capture densely located cells compared with SCCNN. Models for the B-cell panel were developed on a total of 13,545 single-cell annotations collected from 5 of 10 IHC slides and randomly split into training (64%), validation (16%), and testing (20%). Cell annotations of the T-cell and B-cell panels were both generated by A. Akarca and H. Zhang, and reviewed by T. Marafioti. To ensure the quality of cell annotations, we excluded labels disagreed by the two annotators from the dataset.
On average, we identified 69,146 AE 22908 CD8 þ , 109,764 AE 39,028 CD4 þ FOXP3 À and 34,861 AE 28,737 CD4 þ FOXP3 þ in the T-cell panel, achieving an average accuracy of 0.96 on 596 cells. For the B-cell panel, we identified 13,858 AE 13,950 CD20 þ CXCR5 À , 10,399 AE 8,510 CD20 þ CXCR5 þ , and 7,281 AE 4,336 CD79b þ cells with an average accuracy of 0.89 on 2,800 cells. To evaluate the impact of annotation budget on model performance, we retrained the classifier for the B-cell panel with different size of randomly sampled training set, and evaluated the performance on the same hold-out testing set. The model achieved an accuracy of 0.867, 0.879, 0.884, and 0.894 with 25%, 50%, 75%, and 100% of the annotations, indicating that the current amount of annotations was needed to obtain a reliable supervised classifier.

Serial section alignment
T-and B-cell sections in the validation cohort were aligned to the corresponding H&E slide using manually chosen benchmarks. The affine transformation was performed using the "fitgeotrans" function in MATLAB R2019b (RRID:SCR_001622). The quality of registration was evaluated by (i) the Target Registration Error, which represents the distance between original and transformed points normalized by the diagonal length of an image (20); (ii) the concordance between densities of cancer cells and lymphocytes detected on H&E and IHC slides, respectively.
The averaged target registration error of aligned sections was 0.13% and the densities of cancer and lymphocytes detected on IHC significantly correlated to that on H&E sections (cancer: R ¼ 0.69, P ¼ 0.00079; lymphocyte: R ¼ 0.76, P ¼ 8.7 Â 10 À5 ; Supplementary  Fig. S3A and S3B).

Annotation of tertiary lymphoid structures and lymphoid aggregates
Lymphocyte clusters forming a structure were delineated by a senior pathologist on 10 digitized H&E slides from the LATTICe cohort. This manual identification was done while looking at both T-cell and B-cell serial IHC sections (Fig. 4A). Tertiary lymphoid structures (TLS) were marked as the aggregation of lymphocytes in which there was a clear zonality of T-and B-cell components, resulting in a structure. LAGs were marked as lymphocyte clusters without a distinguishable germinal center or admixed with stromal cells. Individual TLS and LAG were assigned to intratumoral, peritumoral, cancer hotspots and non-hotspot regions based on the location of the majority of its region. In total 52 TLSs were identified across 10 slides, with 49 in peritumoral IH, 2 in intratumoral IH, and 1 in non-hotspot regions. 1 out of 52 TLSs had the same area spanning intratumoral and peritumoral IH and was assigned the intratumoral IH. On the other hand, 98 LAGs were identified, with 66 in peritumoral IH, 17 in intratumoral IH, 9 in cancer hotspots and 6 in non-hotspot regions.

Identification of individual immune hotspot
Disconnected IH was identified as cluster of 50Â50 mm 2 grids on the 10 H&E slides of the validation cohort. The spatial connectivity was assessed using the "clump" function in the R package "raster" with directions set as 8. Each individual IH was first segregated and indexed regardless of its distribution, then partially assigned to the peritumoral

Cellular interaction analysis
We calculate the cellular interaction for each individual IH using the coordinates of immune cells within an individual IH as input. The spatially neighboring cells were connected using the Delaunay triangle graph, with the connections exceeding 250 mm excluded. The interactions between each pair of immune cell subsets were defined as the fractions of links between them among all the connections (21).

Data availability
The transcriptomic data and the digital pathology images from the TCGA cohort can be accessed on the GDC portal (portal.gdc.cancer. gov, cohorts TCGA-LUSC and TCGA-LUAD). IHC data related to the validation cohort are available upon reasonable request.

Code availability
Code and data used to generate figures and statistical information can be obtained at https://github.com/yuerua/t_b_spatial_IH. Scripts for implementing the SCCNN pipeline used for cell identification on the H&E and T-cell panel of the validation cohort are available at https://github.com/qalid7/compath. Scripts for the DRDIN pipeline used for cell identification on the B-cell panel of the validation cohort are available at https://github.com/pathdata/UNMaSk/tree/mas ter/HE_cell_detection.

Results
High intratumoral immune hotspot score correlates with poor survival in LUSC Image analysis and spatial statistics were used to identify spatially and histologically distinct immune-rich and cancer cell-rich areas in fresh frozen H&E slides (4)  Quantitative spatial scores S intra/immune , S intra/cancer , and S intra/tissue were used to summarize the amount of intratumoral IH in relation to immune-rich regions, cancer-rich regions, and the entire tissue, respectively (Fig. 1B, see Materials and Methods; ref. 4). We observed a higher S intra/immune in LUAD than LUSC (P ¼ 0.011), and in females than males (P ¼ 0.03). The S intra/cancer was higher in females than males (P ¼ 0.016). No statistically significant difference in the three spatial scores was found for groups categorized by other clinicopathologic characteristics, including age, smoking history, pack years, TMB, and stage (P > 0.05, Table 1).
To determine the prognostic value of spatial scores, we carried out survival analysis in discovery and validation cohorts randomly split from the TCGA LUSC cohort, each containing 231 patients (Supplementary Table S2). Among the three scores, S intra/immune was the only one significantly associated with the overall survival (OS) in both discovery (P ¼ 0.019) and validation (P ¼ 0.043) cohorts ( Fig. 1C; Supplementary Fig. S4A and S4B), with high S intra/immune indicative of an adverse outcome (threshold 0.326 determined in the discovery cohort). To examine the contingency of the prognostic value of S intra/immune on the cohort split, we repeated the splitting of discovery and validation cohort for 100 times. In 36 iterations, S intra/immune showed a significant prognostic effect (P < 0.05) in the validation cohort when dichotomized by the threshold optimized for the discovery cohort, with an averaged P value of 0.224 (SD, 0.261; 95% CI, 0.172-0.275). The percentage of significant iterations was higher than that of other well-known prognosticators including TNM stage (8%), age (4%), pack years (6%), and TMB (20%; Supplementary Table S3; refs. [22][23][24][25], signifying the discriminative ability of S intra/immune to identify patients at high risk. Multivariate analysis showed that S intra/immune remained significant after accounting for age, stage, and smoking pack years in the validation cohort [P ¼ 0.02; HR ¼ 2.0 (1.12-3.7); Fig. 1D]. S intra/immune was also significantly associated with poor outcome when tested separately in patients stratified on the basis of tobacco smoking history available in the TCGA clinical data, differentiating current reformed smokers for less than 15 years, current reformed smokers for more than 15 years, and current smokers (P < 0.05; Supplementary Fig. S5A). However, it was not significant in nonsmokers, potentially due to the small sample size (P ¼ 0.18, n ¼ 16; Supplementary Fig. S5A).
In LUAD tumors, S intra/immune did not show significant prognostic value either in the validation cohort (P ¼ 0.9, n ¼ 236; Supplementary Fig. S6A) or within any of the smoking categories (P > 0.05; Supplementary Fig. S5B). The other two spatial scores S intra/cancer and S intra/tissue were both not significant in the validation cohort (P ¼ 0.094, P ¼ 0.69, n ¼ 236; Supplementary Fig. S6B and S6C).
This differential prognostic value of S intra/immune between LUSC and LUAD pertaining to the spatial heterogeneity of lymphocyte distribution was not related to the overall immune infiltration in both NSCLC subtypes. Lymphocyte abundance, defined as the percentage of lymphocytes detected on the whole slide, was higher in LUAD than LUSC (P ¼ 0.0018). A high level of lymphocyte abundance was associated with favorable OS in LUAD (P ¼ 0.017) but not in LUSC (P ¼ 0.08; Supplementary Fig. S7A and S7B). However, in multivariate analysis adjusted for stage, age and smoking pack years, lymphocyte abundance was not prognostic in either subtype (P > 0.05; Supplementary Fig. S7C and S7D). In addition, S intra/immune did not correlate with lymphocyte abundance in LUSC (R ¼ -0.081, P ¼ 0.083; Fig. 1E). Tumors with high lymphocyte abundance and low S intra/immune displayed accumulation of lymphocytes at the peritumoral region    Supplementary Fig. S8A), whereas in tumors with low lymphocyte abundance and high S intra/immune , the majority of tumor infiltrating lymphocytes were found within the tumor nest ( Supplementary  Fig. S8B). Patients with high lymphocyte abundance but low S intra/immune , indicating a preferential homing of lymphocytes to the peritumoral region, had significantly better prognosis than other patients [P ¼ 0.00035, HR ¼ 0.55 (0.4-0.76), Fig. 1F].
Consistent upregulation of B-cell signatures associated with high intratumoral immune hotspot score To investigate molecular heterogeneity underpinning intratumoral and peritumoral IH, we leveraged RNA-seq data to evaluate differentially expressed genes (DEG) in groups of patients with LUSC stratified by S intra/immune . Gene set enrichment analysis revealed that the Gene Ontology biological processes significantly enriched by these DEGs included immune response, inflammatory response, and activation of immune response ( Fig. 2A). Immune related DEGs were exclusively B-cell markers, including CD79B, MS4A1 (encoding CD20), and CXCR5, which were all upregulated in the high S intra/immune group (Fig. 2B). CXCR5 expression on B cells characterizes the subpopulation migrating into TLS, the organized lymphoid aggregates with morphologically distinct separation of T-cell and B-cell zones and serve as local sites of T-cell priming and B-cell maturation (26)(27)(28)(29)(30). The association between S intra/immune and B-cell infiltration was cross-validated using seven established algorithms to estimate immune cell abundance, including Danaher (11), TIMER (10), CIBERSORT Absolute mode (12), MCP-counter (9), quanTIseq (13), EPIC (14), and xCell (15). Using these methods, we found a compelling pattern for an enrichment of B-cell populations in the high S intra/immune group as compared with the low group according to the median (Fig. 2C).
In contrast with LUSC, S intra/immune high group in LUAD showed a downregulation of B cell-related genes (CD79b, FCER2, CXCR5, MS4A1, CD19; Supplementary Fig. S2B). The enrichment analysis revealed that DEGs between high and low groups in LUAD were enriched for microtubule-based movement and cell-adhesion, unlike the predominance of immune-related pathways in LUSC (Supplementary Fig. S2A).
We next examined the cooccurrence of T-and B-cell genes and signatures in the high S intra/immune LUSC group. B-cell memory genes such as CD27 were significantly correlated with FOXP3 expression (R ¼ 0.8; Supplementary Fig. S9A), compared with a weaker correlation between CD27 and CD8A (R ¼ 0.55; Supplementary Fig. S9B), in line with a previous report (31). Both CD27 and FOXP3 were overexpressed in the high S intra/immune group ( Supplementary Fig. S9A and S9C). The MCP-counter B-lineage signature correlated with the T-cell signature but not the cytotoxic lymphocyte or CD8 T-cell signature, indicating that T-cell subsets more likely to coexist with infiltrating B cells were CD4 þ T cells, and not effector T cells (Supplementary Fig. S9D).
We then investigated the association between S intra/immune and potential indicators of immunotherapy response, including TMB, PD-1, PD-L1 expression, and the neoantigen load (1) in patients with LUSC. We observed a significantly lower TMB in S intra/immune high group (P ¼ 0.032, Supplementary Fig. S10A). On the other hand, there was no statistical difference in PD-1, PD-L1 expression, or number of clonal neoantigens between S intra/immune high and low group (PD-1, P ¼ 0.19; PD-L1, P ¼ 0.12; neoantigen load, P ¼ 0.29; Supplementary  Fig. S10E, S10F and S10G). Consistent with previous results (24), low TMB was associated with adverse outcomes in LUSC (P ¼ 0.014, Supplementary Fig. S10B). Despite the negative association between S intra/immune and TMB, high S intra/immune independently correlated with poor OS in the multivariate analysis including TMB, age, stage, and pack years [P ¼ 0.022, HR ¼ 1.6 (1.07-2.5); Supplementary Fig. S10C]. In addition, high TMB and low S intra/immune codefined a group of patients with significantly prolonged survival (P ¼ 0.0011; Supplementary Fig. S10D).

Abundance of T-and B-cell subpopulations in intratumoral and peritumoral immune hotspots
To further investigate the distribution of B cells in relation to intratumoral and peritumoral IH as well as the interplay between T and B cells, we selected 10 LUSC tumors from the LATTICe cohort, with 6 displaying high S intra/immune (see Materials and Methods). Guided by the gene expression analysis (Fig. 2B), we stained three consecutive sections, respectively, with a panel for B cells (CXCR5/ CD79b/MS4A1, encoding CD20) and neoplastic cells (P40), a T-cell panel (CD4/CD8/FOXP3), and with H&E.
The three sections were spatially aligned using the affine transformation based on manually selected landmarks (averaged misalignment distances: 43.33 mm; averaged target registration error: 0.13%, see Materials and Methods). Deep learning approaches were used to automate the identification of three dominant B-cell subsets (CD20 þ CXCR5 À cells, CD20 þ CXCR5 þ cells, and CD79b þ cells; average accuracy: 0.89) and T-cell subsets (CD4 þ FOXP3 À T cells, CD8 þ T cells, and CD4 þ FOXP3 þ T regulatory (Treg) cells; average accuracy: 0.96, see Materials and Methods). The CD79b þ cells represented B cells with high CD79b expression and different expression levels of CD20 and CXCR5. CD20 þ CXCR5 À and CD20 þ CXCR5 þ cells both showed low expression of CD79b. After the identification of B and T cells, we employed spatial statistics to quantify their distribution and colocalization patterns in intratumoral and peritumoral IH inferred from the corresponding H&E sections (Fig. 3A-E, see Materials and Methods).
In general, there was a higher density of T cells than B cells in the whole tissue, with CD4 þ FOXP3 À and CD8 þ being the first and the second most abundant T-cell subsets ( Supplementary Fig. S3C). There was no significant difference among densities of B-cell subsets (Supplementary Fig. S3C). The majority of T-cell subsets were found in the non-hotspot tissue regions, whereas the B-cell subsets displayed the highest percentages in the peritumoral IH ( Supplementary Fig. S3D). We further compared immune cell compositions between intratumoral and peritumoral IH. Specific B-cell subsets, CD20 þ CXCR5 þ cells and CD79b þ cells, displayed higher densities in peritumoral IH than in intratumoral IH (P < 0.05, Fig. 4B). There was no difference between densities of T-cell subsets in intratumoral and peritumoral IH (P > 0.05, Fig. 4B). The lower B-cell subset densities at the intratumoral IH compared with the peritumoral IH suggested that the concordant increase in B-cell signatures and proportions of intratumoral IH was due to recruitment of B cells to tumors with increased tumorimmune interface, rather than preferential homing of B cells to the intratumoral IH.
TLSs account for a minor proportion of intratumoral and peritumoral immune hotspots CXCR5 and CD79b have been recognized as markers associated with the TLS (28)(29)(30). To investigate the overlap of TLS with intratumoral and peritumoral IH, pathologists manually identified 52 TLSs with prominent germinal centers and 98 lymphoid aggregates (LAG) without a distinguishable germinal center or admixed with stromal cells from 10 H&E slides from the LATTICe cohort (Fig. 4A, see Materials and Methods).
We observed a higher density of CD20 þ CXCR5 þ and CD79b þ cells in TLS within IH than in non-TLS/LAG IH (P < 0.05; Fig. 4E). This is consistent with previous findings of the upregulation of CXCR5 and CD79b in tumors with TLS present (28)(29)(30). TLS also contained a higher amount of CD4 þ T cells, CD4 þ FOXP3 þ Treg cells, and CD20 þ CXCR5 À B-cell subsets compared with non-TLS/LAG IH (P < 0.01, Fig. 4E), whereas only CD20 þ CXCR5 þ and CD79b þ B cells were significantly different between TLS and LAG. Taken together, the majority of IHs are neither TLS nor LAG, based on their different B-and T-cell composition. We henceforth refer to these as non-TLS/LAG IH.

B-and T-cell spatial interactions in non-TLS/LAG immune hotspots
To further investigate the hazard underpinning intratumoral IH while excluding the confounding effect of TLS, we sought to resolve cell spatial interactions that differentiate non-TLS/LAG intratumoral and peritumoral IH. Using the Delaunay triangle graph (21), we developed a method for identifying individual non-TLS/LAG IH based on the connectivity of adjacent regions and construct cellular networks (Fig. 5A). Lymphocyte subset interactions were defined as the proportion of edges connecting two types of lymphocytes (see Materials and Methods). Both interactions and compositions of lymphocyte subsets were significantly more diverse in the peritumoral IH than in the intratumoral IH, as indicated by the higher Shannon index (P ¼ 0.003; P ¼ 0.003; Fig. 5B). Compared with intratumoral IH, peritumoral IH had higher frequencies of interactions between CD20 þ CXCR5 þ and CD4 þ FOXP3 þ , CD4 þ FOXP3 À , but not CD8 þ cells (univariate logistic regression P < 0.05; Fig. 5C). On the other hand, none of the interactions involving CD79b þ B cells were observed to be significant despite the enrichment of CD79b þ B cells in peritumoral IH (P > 0.05, Fig. 5C), suggesting an independence between cell interactions and abundance.
High abundance of B cells in non-TLS/LAG intratumoral immune hotspots was linked to a decrease in the CD8 þ /CD4 þ FOXP3 þ ratio We assessed the localized T-cell regulation by the CD8 þ / CD4 þ FOXP3 þ ratio within individual non-TLS/LAG IH (see Materials and Methods). There was no difference between the CD8 þ /CD4 þ FOXP3 þ ratio at intratumoral and peritumoral IH (intratumoral IH: mean ¼ 4.66, SD ¼ 7.09, peritumoral IH: mean ¼ 3.72, SD ¼ 4.3, P ¼ 0.53). To further dissect the role of B cells on T-cell regulation within non-TLS/LAG IH, each non-TLS/LAG intratumoral and peritumoral IH was classified on the basis of its CD8 þ /CD4 þ FOXP3 þ ratio. IH with low CD8 þ /CD4 þ FOXP3 þ ratio had significantly higher percentages of CD20 þ CXCR5 À and CD20 þ CXCR5 þ B cells (defined as the number of cells divided by the total count of 6 immune cell types in each IH) compared to IH with high CD8 þ /CD4 þ FOXP3 þ ratio, regardless of their location (P < 0.05, Fig. 6A).
We then investigated whether the reduced CD8 þ /CD4 þ FOXP3 þ ratio at B-cell-enriched IH was associated with the exclusion of CD8 þ T cells or the recruitment of Treg cells. B-cell-enriched intratumoral and peritumoral IH had significantly lower percentages of CD8 þ T cells compared with IH with low B-cell percentages (intratumoral IH: Fig. 6B). In contrast, CD4 þ FOXP3 þ T-cell percentage was not correlated with B-cell abundance in intratumoral IH or peritumoral IH (P > 0.5, Fig. 6B). In addition, we observed a strong negative correlation between CD8 þ and CD4 þ FOXP3 þ percentages in intratumoral (P ¼ 2.8 Â 10 À5 ) but not in peritumoral IH (P ¼ 0.060, Fig. 6C). Taken together, these findings suggest that while a generic pattern of decreased CD8 þ cells but not Treg cells in B-cell-enriched IH is observed regardless of their location, Treg cells may have a bigger impact on the regulation of CD8 þ cells in intratumoral IH compared with peritumoral IH, due to the presence of negative correlation between CD8 þ T-cell and Treg abundance exclusively in intratumoral IH.

Discussion
Heterogeneous spatial organization of tumor-infiltrating lymphocytes was associated with immune evasion and tumor evolution in NSCLC (2,33). Specifically, the presence of TLS has been shown to correlate with prolonged survival in NSCLC and a variety of cancer types (28,(33)(34)(35)(36). However, immune rich areas beyond TLS are not well studied. Immune-rich areas differ from TLS in terms of cellular densities, organizations, and spatial distributions (37)(38)(39). Understanding the role of these immune rich areas in cancer progression and treatment response may be crucial in improving therapies involving immune checkpoint blockade and induction of TLS against cancer (40). In this study, we integrated transcriptional profiling and computational pathology to deconvolve the spatial heterogeneity of immune rich areas identified by the hotspot analysis in NSCLC. We measured the proportion of IH colocalizing with tumor as a score S intra/immune . The score was a prognostic factor associated with unfavorable outcomes in LUSC independent of age, stage, and pack-years. Despite S intra/immune being higher in LUAD than LUSC, the score did not show prognostic value or association with immune-related genes in LUAD. Such discrepancy adds to the current paradigm of the distinct biology and immune evasion mechanisms of the two histologic subtypes, such as the different prevalence of HLA LOH and clonal neoantigens (23). This finding highlights the importance of studying IHs beyond TLS and indicates that the prognostic value of lymphocyte clustering is nonexclusive to breast cancer.
The immune spatial score S intra/immune may be a B-cell-enriched spatial signature, given the consistent upregulation of B-cell signatures in the TCGA cohort. It was later confirmed in the IHC experiments in the validation cohort based on B-cell accumulation at the IHs. In addition, LUSC with high S intra/immune displayed a lower TMB, raising the possibility that increased IH located inside the tumor nest might characterize tumors undergoing clonal selection or adaptation to the immune microenvironment. One of the potential reasons for S intra/immune being an adverse survival indicator is that S intra/immune -high tumors had reduced peritumoral IH and therefore low TLS abundance, which has been associated with adverse outcome in NSCLC (33). TLS are characterized as organized lymphoid aggregates with a clear separation of T-and B-cell zones (28,29). In the validation cohort, we found TLS frequently located outside the tumor nest, and contain higher densities of B cells and CD4 þ T cells than intra and peritumoral IH. The majority of peritumoral and intratumoral IH were immune-rich areas without organized lymphoid structure. Given the minor proportion of TLS within both intratumoral and peritumoral IH, it is unlikely that S intra/immune is a mere reverse approximate of TLS abundance. Comparisons between non-TLS/LAG intratumoral and peritumoral IH showed increased interactions between CD20 þ CXCR5 þ B and CD4 þ T cells, Treg cells in the non-TLS/LAG peritumoral IH. B cells have been shown to present antigen and to activate CD4 þ FOXP3 À T cells at the interfollicular region in lymph nodes and within TLS (41,42). It is possible that CD20 þ CXCR5 þ B cells in non-TLS/LAG peritumoral IH also exhibit such proinflammatory functions, contributing to the antitumor immune response in peritumoral IH. The enhanced interactions between CD20 þ CXCR5 þ B cells and Treg cells suggested that B cells in peritumoral IH were subjected to stronger regulation, potentially due to the accumulation of tumorspecific autoantibodies in peritumoral IH (43). Furthermore, both non-TLS/LAG peritumoral and intratumoral IH exhibit a decreased CD8 þ T-cell percentage and CD8 þ /CD4 þ FOXP3 þ ratio in response to the increased proportions of CD20 þ CXCR5 À and CD20 þ CXCR5 þ B cells. In addition, an increase in Treg percentage was associated with decreased CD8 þ T-cell proportions in non-TLS/ LAG intratumoral IH but not in peritumoral IH (Fig. 6C). This negative relationship in intratumoral IH may indicate an immune escape strategy of tumor cells, where some tumor islands favor the recruitment of Treg over cytotoxic T cells, leading to a low-CD8 þ high-Treg intratumoral IH phenotype (44,45). Taken together, these results suggest an immunosuppressive TME imposed by the accumulation of B cells and T-cell regulation, with cytotoxic T cells subjected to greater regulation in intratumoral IH than in peritumoral IH. Therapeutic agents targeting B cells and Treg cells, such as rituximab and daclizumab, might skew the balance between intratumoral and peritumoral IH toward an immunogenic microenvironment.
This study has several limitations. The IHC panel was designed based on the DEGs inferred from the TCGA cohort, which did not include markers for the functionality of T cell subsets such as PD-1, ICOS, and GZMB (37). It is possible that CD8 T cells residing in intratumoral IH are prone to exhaustion and are expected targets of immune checkpoint blockade therapy. Likewise, this panel was unable to reveal the potential immunosuppressive role of B cells. For example,  in vitro experiments demonstrated that B regulatory cells promoted the expression of FOXP3 and CTLA4 in Treg via secretion of TGFb and cell-cell contact (46). In addition, na€ ve B cells and follicular lymphoma B cells have been shown to induce the conversion of CD4 þ T cells to Treg in a cell-cell contact manner (47,48), which potentially contributed to the increased signal of Treg observed in the TCGA transcriptomic data. Also, B cells in intratumoral IH can possibly acquire an IL10-secreting phenotype that displays an immunoregulatory role (49). To fully appreciate the components, functionalities, and spatial interactions of lymphocytes and myeloid cells in intra and peritumoral IH, future studies are warranted. Also, the heterogeneous distribution of IH with respect to tumor glands might be attributed to a diverse immune evasion capability of tumor clones, as indicated by the branching evolution and tissue architecture of invasive glandular tumors (50). Another caveat is the identification of TLS, which was dependent on the variable choices of markers and morphologic criteria across studies (33). Also, the shape and structure of TLS may be varied by the sectioning angle. As a result, despite best efforts, the annotated TLS and LAG regions in the validation cohort may exclude partial TLS regions and immature TLS. However, the presence of CD20 þ CXCR5 þ and CD79b þ was consistently observed in non-TLS/LAG IH, indicating that the two types of B cells preferentially reside in the peritumoral IH with respect to intratumoral IH regardless of the presence of TLS. In addition, further investigation on the distribution of IHs in preinvasive lesions has the potential to benefit the patient risk assessment, as immune surveillance in carcinoma in situ has been associated with lesion regression (51).
In conclusion, we propose a new spatial signature, S intra/immune , of protumor response orchestrated by B cells and T cells, as the balance between intratumoral and peritumoral IH (Fig. 7). Our study demonstrated a negative association between S intra/immune and OS in LUSC. Increased S intra/immune was correlated with an upregulated B-cell signature. Although limited by sample size, our preliminary results indicated the accumulation of CD20 þ CXCR5 À /CD20 þ CXCR5 þ / CD79b þ B cells and the minor proportion of TLSs in IH. The interplay between CD8 þ T cells, Tregs, and B cells underlying the spatial clustering pattern of lymphocytes characterized a protumor micro-environment in the intratumoral IH that may fuel the progression of LUSC.
Authors' Disclosures D.A. Moore reports personal fees from AstraZeneca, Eli Lilly, Takeda, ThermoFisher, Amgen, Janssen, MIM Software, and Bristol-Myers Squibb outside the submitted work. K.S.S. Enfield reports other support from European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement no. 838540 and grants from The Royal Society (RF\ERE\210216) during the conduct of the study, and grants from Bristol Myers Squibb outside the submitted work. M. Jamal-Hanjani reports grants from Cancer Research UK during the conduct of the study; grants from Cancer Research UK, Rosetrees Trust, NIHR, NIH National Cancer Institute, IASLC International Lung Cancer Foundation, and Lung Cancer Research Foundation; nonfinancial support from Achilles Therapeutics Scientific Advisory Board and Steering Committee; other support from Pfizer, Astex Pharmaceuticals, Oslo Cancer Cluster, and Bristol Myers Squibb outside the submitted work. R. Salgado reports other support from BMS, Exact Sciences, Roche, Merck, Puma, and BMS outside the submitted work. N. McGranahan reports a patent for PCT/ EP2016/059401 issued, a patent for PCT/EP2016/071471 issued, and a patent for PCT/GB2018/052004 issued. C. Swanton reports grants and personal fees from AstraZeneca, Boehringer-Ingelheim, Bristol Myers Squibb, Invitae (previously Archer DX), Ono Pharmaceuitcal, Personalis, Medicxi, China Innovation Centre of Roche (CICoR), Sarah Cannon Research Institute, Metabomed, Amgen, Novartis, GSK, MSD, Illumina, and Roche-Ventana; other support from Bicycle Therapeutics, Apogen Biotech, and personal fees and other support from Epic Bioscience during the conduct of the study. In addition, C. Swanton has a patent for PCT/2017/028013 for methods to lung cancer; a patent for PCT/EP2016/ 059401 targeting neoantigens pending, a patent for identifying patent response to immune checkpoint blockade (PCT/EP2016/071471) pending, a patent for determining HLA LOH (PCT/GB2018/052004) pending, a patent for predicting survival rates of patients with cancer (PCT/GB2020/050221) pending, a patent for identifying patients who respond to cancer treatment (PCT/GB2018/051912) pending, a patent for methods for lung cancer detection (US20190106751A1) pending, and a patent for (PCT/GB2017/053289) relating to assay technology to detect tumour recurrence pending. C. Swanton is an AstraZeneca advisory board member and chief investigator for the AZ MeRmaiD 1 and 2 clinical trials and is also co-chief investigator of the NHS Galleri trial funded by GRAIL and a paid member of GRAIL's Scientific Advisory Board. C. Swanton is a co-founder of and has stock options in Achillles Therapeutics. Y. Yuan reports personal fees from Roche and Merck and Co Inc. during the conduct of the study. No disclosures were reported by the other authors.

Cancer hotspots
Peritumoral immune hotspots Intratumoral immune hotspots Figure 7. Immune spatial balance in LUSC. Spatial positioning of immune rich areas can reflect an immunosuppressive or an immunostimulant microenvironment. Increased intratumoral IH compared with peritumoral IH may reflect an immunosuppressive TME, characterized by decreased TLS abundance, lower densities of CD20 þ CXCR5 þ , and CD79b þ B-cell densities, decreased interactions between CD20 þ CXCR5 þ B cells and CD4 þ T cells, and the presence of a negative correlation between CD8 þ T-cell and Treg cell abundance. Also, high ratio of intratumoral IH is coupled with increased B-cell and Treg infiltration in the whole tumor. Cancer hotspots, peritumoral IHs, and intratumoral IHs are denoted by red, blue, and green colors, respectively. An IH can consist of peritumoral and intratumoral compartments, denoted as circles with a combination of green and blue.