Geospatial immune variability illuminates differential evolution of lung adenocarcinoma

Abstract

Remarkable progress in molecular analyses has improved our understanding of the evolution of cancer cells toward immune escape1-5. However, the spatial configurations of immune and stromal cells, which may shed light on the evolution of immune escape across tumor geographical locations, remain unaddressed. We integrated multiregion exome and RNA-sequencing (RNA-seq) data with spatial histology mapped by deep learning in 100 patients with non-small cell lung cancer from the TRACERx cohort6. Cancer subclones derived from immune cold regions were more closely related in mutation space, diversifying more recently than subclones from immune hot regions. In TRACERx and in an independent multisample cohort of 970 patients with lung adenocarcinoma, tumors with more than one immune cold region had a higher risk of relapse, independently of tumor size, stage and number of samples per patient. In lung adenocarcinoma, but not lung squamous cell carcinoma, geometrical irregularity and complexity of the cancer-stromal cell interface significantly increased in tumor regions without disruption of antigen presentation. Decreased lymphocyte accumulation in adjacent stroma was observed in tumors with low clonal neoantigen burden. Collectively, immune geospatial variability elucidates tumor ecological constraints that may shape the emergence of immune-evading subclones and aggressive clinical phenotypes.

Trial registration: ClinicalTrials.gov NCT01888601.

Conflict of interest statement

Competing Interests

Y.Y. has received speakers bureau honoraria from Roche and is a consultant for Merck and Co Inc. C.S. receives grant support from Pfizer, AstraZeneca, BMS, Roche-Ventana, Boehringer-Ingelheim and Ono Pharmaceutical. C.S. has consulted for Pfizer, Novartis, GlaxoSmithKline, MSD, BMS, Celgene, AstraZeneca, Illumina, Genentech, Roche-Ventana, GRAIL, Medicxi, and the Sarah Cannon Research Institute. C.S. is a shareholder of Apogen Biotechnologies, Epic Bioscience, GRAIL, and has stock options in and is co-founder of Achilles Therapeutics. M.A.B. is a consultant for Achilles Therapeutics. S.L. receives research funding to her institution from Novartis, Bristol Meyers Squibb, Merck, Roche-Genentech, Puma Biotechnology, Pfizer, Eli Lilly and Seattle Genetics. S.L. has acted as consultant (not compensated) to Seattle Genetics, Pfizer, Novartis, BMS, Merck, AstraZeneca and RocheGenentech. S.L. has acted as consultant (paid to her institution) to Aduro Biotech, Novartis, and G1 Therapeutics. D.A.M. has received speaker’s fees from AstraZeneca. M.J.H. is a member of the Advisory Board for Achilles Therapeutics.

Figures

Extended Data Fig. 1. CONSORT diagrams for TRACERx 100 and LATTICe-A histology cohorts and patient characteristics.

a. TRACERx CONSORT diagram to illustrate sample collection and analysis of regional and diagnostic histology samples, as well as the overlap with RNA and DNA studies. b. TRACERx patient characteristics for the histology cohort. c. LATTICe-A CONSORT diagram (= 970 LUAD patients). Legends for ‘type of the analysis’ correspond to panel a. d. Demographics and clinical patient characteristics for TRACERx (top three panels) and LATTICe-A (bottom three panels) showing the distribution of age (colored by sex), distribution of smoking pack years and the proportion of patients in each pathological stage. Horizontal lines indicate the median value.

Extended Data Fig. 2. Validation of the automated single-cell classification for H&E.

a. A scatter plot showing the correlation between H&E-based adjacent-tumor lymphocytes/stromal and pathology TIL estimates in diagnostic samples (= 98 diagnostic slides/patients). b. Scatter plots showing the correlations between H&E-based tumor cellularity estimate and ASCAT/VAF purity scores (= 238 regions; 83 patients). c. A scatter plot showing the correlation between H&E-based estimate of lymphocyte percentage among all cells and RNA-seq-based CD8+ signature using the Danaher et al. method (= 142 regions; 56 patients). d. A scatter plot showing the correlation between H&E-based estimate of lymphocyte percentage among all cells and CD8+ cell percentage in IHC in the diagnostic samples (= 100 diagnostic slide/patients). e. Scatter plots showing the correlation between H&E-based lymphocyte percentage versus pathological scores of overall lymphocytic cell fraction, and adjacent-tumor lymphocytes/stromal versus pathology TIL estimates in an external cohort (LATTICe-A, = 80 diagnostic slides/patients). f. Illustrative example to show the spatial alignment of TTFl/CD45/SMA-stained IHC and H&E images obtained using sequential staining on the same tissue microarray section for biological validation, g. A scatter plot showing the correlation between stromal cell percentage determined by H&E and SMA+ cell percentage per LUAD image tiles of size 100μm2 (= 144). The experiment was conducted once using one TMA (= 33 cores/patients). The shading indicates 95% confidence interval.

Extended Data Fig. 3. Distribution of regional lymphocytic infiltration according to pathological stage.

All available patients’ data have been used in this figure except for the standard deviation tests excluding patients with a single tumor region. Patients without pathological staging information from the LATTICe-A cohort were also removed. a, b, c, top row: TRACERx and bottom row: LATTICe-A. Horizontal lines indicate the median value. a. Distribution of the standard deviation of regional lymphocyte percentage for LUAD and LUSC patients in TRACERx (= 69), and LUAD in LATTICe-A (= 814). b. Distribution of the standard deviation of regional lymphocyte percentage across pathological stages (= 69 for TRACERx, 814 for LATTICe-A). c. Distribution of regional mean of lymphocyte percentage across stages (= 79 for TRACERx, 827 for LATTICe-A). d. No significant difference among stages with respect to standard deviation (= 69 for TRACERx, 814 for LATTICe-A) or mean (= 79 for TRACERx, 827 for LATTICe-A) of regional lymphocytic infiltration. Left panel, TRACERx and right panel, LATTICe-A. Correction for multiple testing was applied in d, for each cohort individually. A two-sided, non-parametric, unpaired, Wilcoxon signed-rank test was used; each dot represents a patient; the mean value is annotated with a large dot; the median value is represented by a thick horizontal line; minimum and maximum values are indicated by the extreme points; the first and third quantiles are represented by the box edges; and the violin shape shows the data distribution as a kernel density estimation.

Extended Data Fig. 4. Validation of immune phenotype classification.

a. The proposed immune classification imposed on density plot showing distribution of lymphocyte percentage. The middle zone corresponds to the intermediate phenotype, red zone for immune hot and blue zone for immune cold. Black dash line shows the median. This classification was validated after applying small perturbations to the thresholds to reclassify regional immune phenotypes, illustrated as grey dash lines: no intermediate zone (i.e. hard median for separating hot and cold), standard deviation (SD)/2 above and below the median, SD/3 and SD/6. b. Forest plots to show repeated multivariate Cox regression tests for the number of immune cold regions using these new classifications (= 79 patients), after accounting for stage, total number of samples, upper quartile of clonal neoantigens determined for LUAD and LUSC individually, and other clinical parameters. Box plots showing difference in genomic distance for pairs of hot regions compared with pairs of cold regions for LUAD and LUSC separately (LUAD: = 45 hot pairs, 45 cold pairs for no intermediate zone; = 19 hot, 25 cold for SD/2; = 25 hot, 33 cold for SD/3; = 32 hot, 41 cold for SD/6. LUSC: = 32 hot pairs, 54 cold pairs for no intermediate zone; = 19 hot, 27 cold for SD/2; = 19 hot, 37 cold for SD/3; = 27 hot, 41 cold for SD/6.). c. Box plots showing significant difference in CD8+ RNA-seq signature using the Danaher method between regions of hot and cold phenotype across all classification schemes (= 219 for SD/4; 275 for no intermediate zone; 173 for SD/2; 204 for SD/3; 237 for SD/6). d. Distribution and difference of lymphocytic infiltration for LUAD versus LUSC regions in TRACERx (= 275 regions; 85 patients) as well as distribution for LUAD in LATTICe-A (= 4,324 samples; 970 patients). Horizontal lines in the distribution plots indicate mean values. For statistical comparisons among groups, a two-sided, non-parametric, unpaired, Wilcoxon signed-rank test was used, unless stated otherwise.

Extended Data Fig. 5. Concordance between histology deep learning and RNA-seq immune classification.

a. A box plot showing the difference in pathology TIL estimates between immune hot and immune cold regions (= 219). Pathology TIL estimates score fraction of stroma containing TILs, whereas immune classification was defined based on the percentage of lymphocytes in all cells within a slide, b. A confusion matrix to compare RNA-seq and deep learning histology immune classifications (discarding immune intermediate regions, = 109 regions (57 LUAD, 37 LUSC, 15 other histology subtypes); 52 patients). The p-value was generated using a two-sided Fisher’s exact test for overlap. c. A box plot showing the difference in the fraction of immune hotspots in regions where the two classifications are in agreement (= 78; labeled as ‘In agreement’) against the discrepant regions (= 31, labeled as ‘Discrepant’). Each dot represents a region, the median value is indicated by a thick horizontal line; minimum and maximum values are indicated by the extreme points; and the first and third quantiles are represented by the box edges. d. Box plots to support the overall consistency between H&E-deep learning and RNA-seq methods by comparing different immune scores as well as ASCAT tumor purity between immune hot/high and cold/low tumor regions (all R-values < 0.0001). Top row, H&E-deep learning immune classification (= 219; except the ASCAT purity box plot = 186 regions), bottom row, RNAseq derived immune classification (= 142; except the ASCAT purity box plot, = 141 regions). For statistical comparisons among groups, a two-sided, non-parametric, unpaired, Wilcoxon signed-rank test was used, unless stated otherwise.

Extended Data Fig. 6. Genomic and survival analysis of tumor regions according to immune phenotypes.

a. A box plot showing the difference in genomic distances for pairs of immune hot versus immune cold regions within the same LUSC patients (= 59 pairs). A two-sided, non-parametric, unpaired, Wilcoxon signed-rank test was used. b. Forest plots to show the univariate prognostic value for the number of immune low regions (both as continuous and dichotomized at the median (≤1 versus >1)), or the number of immune high regions, using the immune classification generated by RNA-seq-based infiltrating immune cell populations in 64 TRACERx tumors (41 LUAD, 16 LUSC and 7 other histology subtypes). c. Forest plots showing multivariate Cox regression analyses in both TRACERx (= 79 patients; LUAD and LUSC combined) and LATTICe-A (= 651 LUAD patients representing a subset with complete stage and smoking pack years data) with the number of immune cold regions dichotomized at the median (≤1 versus >1). This remains significant when the number of immune cold regions was replaced as a continuous variable, in the same multivariate model, (R = 0.019 in TRACERx and < 0.001 in LATTICe-A, for the number of immune cold regions). Clonal neoantigens were dichotomized using the upper quartile, determined individually for LUAD and LUSC tumors. d. The same test in c when tumor size (in mm) was also controlled in the multivariate model in LATTICe-A. This test also remained significant for a bigger group of patients with complete stage data, but missing pack years information (= 815, R < 0.001, HR = 1.4[1.1-1.8]). e. Forest plots to compare the prognostic value of regional immune scores as well as diagnostic H&E and IHC scores for relapse-free survival in TRACERx (= 79 patients, LUAD and LUSC combined). Wherever possible, these immune features were tested in LATTICe-A (= 970 patients). To compare the prognostic value of the number of immune cold region with other immune features, LATTICe-A comparisons were conducted in Cox multivariate regression models to include every immune feature after correcting for the number of immune cold regions in the same model. Each variable’s HR is plotted with a 95% confidence interval; all R-values were adjusted for multiple testing; and the size of the circles denotes–log10(R). For the sake of visualization, a minor adjustment was made to the HR for the number of cold regions/total number of regions in LATTICe-A from 0.88[0.57-1.3] to 0.99[0.97-1.3], SD: standard deviation, used for measuring variability of lymphocyte percentage among samples within a tumor, f. Forest plots using Cox multivariate regression analysis showing that the prognostic value of the number of immune cold regions was independent of: 1) genetic measure, subclonal copy number alteration (obtained from ); 2) tumor cellularity from DNA-seq-based ASCAT purity, 3) tumor cellularity measured by deep learning-based cancer cell percentage, g. Kaplan Meier curves to illustrate the difference in relapse-free survival for TRACERx patients including other histology types (= 85; representing all TRACERx patients in the multiregion histology cohort) with high and low number of immune cold regions, dichotomized by its median value. Log-rank R = 0.0017. h. Forest plot using Cox regression for the multivariate survival analysis for the number of immune cold regions in TRACERx including patients with other histology subtypes (= 85).

Extended Data Fig. 7. Fractal dimension and relationships with stromal cells.

a. Distribution of the average minimum Euclidean distance between a stromal cell to its neighboring cancer cell. For every stromal cell in a tumor region slide, the minimum distance to nearest cancer cell was computed. This distance was then averaged for all identified stromal cells in every region to plot the distribution (= 275 regions; 85 patients). b. Distribution of the fractal dimension of the cancer-stroma cell interface for histology types in the TRACERx cohort (= 275 regions; 85 patients). c. Box plots to show the difference in fractal dimension between immune hot and cold regions in TRACERx LUAD (= 113) and LUSC (= 84). d. Box plots showing the difference in stromal cell percentage between immune hot and cold regions in all (= 219), LUAD (= 113), and LUSC (= 84). e. Scatter plots showing the correlation between fractal dimension and percentage of cells that are stromal or cancer in all tumor regions (= 275 regions; 85 patients). This shows that fractal dimension was independent of tumor cell composition, with only a weak correlation with stromal cell percentage and no correlation with tumor cellularity. f. Box plots showing the difference in fractal dimension between LUAD tumor regions harboring an LOH event for HLA type A (= 106), type B (= 113), type C (= 108) versus regions that do not, adjusted for multiple comparisons with the corresponding test in Fig. 4c. g. The same test in f repeated for LUSC tumor regions (= 87) for HLA of any type. h. Box plots showing the difference in tumor-level fractal dimension using the maximum value of regional measures between LUAD tumors (= 48) harboring a single LOH event for any HLA type, HLA type A, type B and type C versus tumors that do not, independent of predicted clonal neoantigens. Each p-value was generated using a multiple regression linear model and was also adjusted for multiple testing correction. i. The same test in h repeated for LUSC tumors (= 29) for HLA of any type. For statistical comparisons among groups, a two-sided, non-parametric, unpaired, Wilcoxon signed-rank test was used, unless stated otherwise.

Extended Data Fig. 8. Relationship of immune subsets and spatial TILs in LUAD.

a. Spearman’s correlations between immune scores in diagnostic slides and genetic measures including predicted neoantigens and HLALOH in LUAD patients (= 46). ITLR: intra-tumor lymphocytes to total tumor cell ratio. Only significant correlations after multiple testing are highlighted (rho = 0.37, R = 0.035). b. Examples of registered H&E and IHC tiles. The green cross denotes a manually placed landmark repeated 238 times on pairs of H&E-IHC image tiles. The Euclidean distance (difference in, coordinates) was computed between the two landmarks which was then c. shown as a distribution to represent the accuracy of the registration (= 249 total H&E-IHC image tiles, maximum five landmarks per a pair of tiles). The average distance between matching landmarks was 9.57μm and the distribution is within the expected range of maximum distance between four serial sections (16μm). d. Box plots to illustrate the difference in percentage of immune cell subsets among adjacent, intra and distal-tumor lymphocytes (= 20 image tiles), a non-parametric, paired Wilcoxon test was used.

Extended Data Fig. 9. Summary of immune and genomics features in NSCLC.

An extended heatmap showing all immune variables described in TRACERx across all patients (= 275 regions; 85 patients), along with genetic measures and clinical parameters. Each column represents a tumor, grouped by their histologic subtype. Tumor regions (illustrated as dots) were assigned to immune hot, immune cold and intermediate phenotypes based on percentage of lymphocytes in all cells following H&E-based deep learning analysis. Cancerstromal fractal dimension, defined using the maximum fractal dimension in regions of a patient, using the median as cut-off to determine high and low groups.

Figure 1. The computational pathology deep learning pipeline for dissecting heterogeneous NSCLC tumor microenvironment.

a. Histology sample generation in Lung TRACERx. To preserve morphology and generate good quality histology, samples from the same tumor regional frozen blocks specifically collected for TRACERx and generated molecular data, were re-embedded in formalin fixed paraffin (FFPE). From these, H&Estained tumor section slides were generated. In addition, H&E section and triplex CD4/CD8/FOXP3 IHC slides were also generated from diagnostic blocks that represent clinical standard sampling, b. Our multistage deep learning pipeline consists of three key stages: fully automated tissue segmentation, single-cell detection and classification. The final output is shown as an image with all cells identified. For more details, please see the ‘Training the deep learning pipeline’ section of the Methods. c. Illustrative 3-dimensional distribution of input image patches in the feature space learned by the convolutional neural networks, using Principal Component Analysis. The feature clusters were pseudocolored to display segregation for four cell types in H&E, and d CD8+, CD4+FOXP3+, CD4+FOXP3- and “other” cell class (hematoxylin cells) in IHC, respectively, e. The deep learning single-cell classification model was trained using expert pathology annotations from a variety of TRACERx samples (diagnostic, regional, TMA). The trained model was then applied to the remaining TRACERx samples (predominantly LUAD and LUSC) and the LATTICe-A cohort (only LUAD), identifying over 171 million cells in TRACERx and over 4.9 billion cells in LATTICe-A. WSI: whole-section image, f. Biological validation of the deep learning approach. H&E and IHC images generated from the same TMA slide were virtually integrated for comparison of H&E-based cell classification and cell type marker expression. For each marker, the experiment was conducted once using a single TMA (cores/patients = 48 TTF1; 38 CD45). Scale bars represent 100μm. g-h. Correlations between cancer/lymphocyte cell percentage determined by H&E and TTF1+ (tumor marker)/CD45+ (immune marker)cellpercentage per LUAD image tilesof size 100μm2 ( = 100 TTF1; 83 CD45). The shading indicates 95% confidence interval.

Figure 2. Geospatial heterogeneity of lymphocytic infiltration in the TRACERx cohort.

a. Representative examples of immune hot and immune cold multi-region H&E samples, scale bars represent 100μm. b. Each column represents a tumor, grouped by their histologic subtype (the “Other” group consists of adenosquamous carcinoma, large cell neuroendocrine carcinoma, pleomorphic carcinoma, and sarcomatoid carcinoma of pleomorphic type arising from adenocarcinoma). Tumor regions (illustrated as dots) were assigned to immune hot, immune cold, and intermediate phenotypes based on percentage of lymphocytes in all cells following H&E-based deep learning analysis. CD8+/CD4+FOXP3/CD4+FOXP3+ percentages based on automated analysis of the IHC diagnostic samples are also shown. c. A heatmap showing gene expression patterns of 14 immune cell populations across tumor regions, each row represents a tumor region (= 142). The three clusters correspond to the proposed immune regional classification as shown in b. d. Significant enrichment of all immune cell populations in hot regions, as compared to cold regions, particularly for the immune activating cell subsets, including cytotoxic, B-cell, and natural killer cells (= 109 regions; 52 patients). A two-sided, nonparametric, unpaired, Wilcoxon signed-rank test was used for each box plot, all R-values were corrected for multiple comparisons. Thick horizontal lines indicate the median value; outliers are indicated by the extreme points; the first and third quantiles are represented by the box edges; and vertical lines indicate the error range.

Figure 3. Evolution of immune escape, and survival analysis in TRACERx and LATTICe-A.

a. A box plot showing the difference in genomic distances for pairs of immune hot or immune cold regions within the same patients in LUAD (= 66 pairs), b. A box plot showing the difference in mutational distance between the dominant subclones in pairs of immune hot or immune cold regions via their last common ancestor in LUAD (= 23 immune cold pairs; 15 immune hot pairs). This distance was calculated by taking the furthest dominant clone (cancer cell fraction (CCF) ≥ 75%) from the trunk, and it remained significant when the dominant clone closest to the most recent common ancestor of each tree was considered (R = 0.02). c. Illustrative examples of tumor phylogenetic trees for a pair of immune hot and immune cold regions. Dominant subclones were labelled and their last common ancestor (annotated with arrows) was then identified. Minor (CCF < 75%) or undetected clones were neglected in this analysis. d,e. Kaplan-Meier curves illustrating the difference in disease-free survival according to the number of immune cold regions, dichotomized by the median value, in TRACERx (d) (LUAD and LUSC, = 79 patients, 249 regions) and LATTICe-A (e) (LUAD, = 970 patients, 4,324 samples). The same deep learning histology analysis and immune regional classification developed for TRACERx were applied directly to LATTICe-A. WSI: whole-section image, f. Forest plots showing multivariate Cox regression analyses in TRACERx (= 79 patients; LUAD and LUSC). Clonal neoantigens were dichotomized using the upper quartile, determined individually for LUAD and LUSC tumors. g. Forest plots showing multivariate Cox regression analyses in LATTICe-A (= 651 LUAD patients with complete stage and smoking pack years data). For the patient subset with complete stage data but missing pack years information, the test remained significant (= 827, R < 0.001, HR = 1.4[1.1-1.9]). For statistical comparisons among groups, a two-sided, non-parametric, unpaired, Wilcoxon signed-rank test was used, unless stated otherwise.

Figure 4. Association of spatial histology with genetic alterations relevant to immune surveillance.

a. An illustrative example of fractal dimension calculated by the box-counting algorithm to quantify the geospatial complexity of the cancer cell-stromal cell interface. By examining boxes of decreasing sizes that contain both cancer and stromal cells, the box counting algorithm quantifies the rate at which the geometrical details of cancer-stromal interface develop at increasingly fine scales. Blue box illustrates the smallest box of 20μm by 20μm in size. Scale bar represent 100μm. An example of a fractal structure displaying geometrical self-similarity is shown below the panel, b. A box plot to illustrate the significant difference in fractal dimension between all TRACERx immune hot and cold regions (= 219). c. A box plot showing a significant difference in fractal dimension between LUAD tumor regions (= 116) harboring an LOH event for class 1 HLA of any type versus regions that do not, adjusted for multiple comparisons with the remaining HLA type-specific tests (see Extended Data Fig. 7f). d. Illustration of the adjacent-tumor lymphocyte/stroma ratio inferred by spatial modeling of cancer cell density (contours) and lymphocyte classification into spatial compartments. Cell classification in IHC sample of the same block was shown for comparison. Scale bars represent 50μm. e. A box plot showing the difference in the adjacent-tumor lymphocyte/stroma ratio between high (≥ median) and low (< median) clonal neoantigens for all LUAD patients in TRACERx (= 61). f. Illustration of image registration to spatially align serial sections of H&E and IHC and generate a virtual composite map of T cell subset in the context of cancer/stroma density. T cell subsets classified in the IHC were projected onto the cancer density map inferred from H&E, so that they can be classified into adjacent-tumor, intra-tumor, and distal-tumor compartments, g. A box plot showing significantly higher ratio of CD8+ to CD4+FOXP3+ cells in adjacent-tumor and intratumor lymphocytes compared with distal-tumor lymphocytes in registered LUAD image tiles (= 20 image tiles, using paired Wilcoxon test). For statistical comparisons among groups, a two-sided, non-parametric, unpaired, Wilcoxon signed-rank test was used, unless stated otherwise.