Integrating AI-Powered Digital Pathology and Imaging Mass Cytometry Identifies Key Classifiers of Tumor Cells, Stroma, and Immune Cells in Non–Small Cell Lung Cancer

The cancer diagnostic workflow is commonly performed through the visual examination by a trained pathologist of histopathologic slides stained with hematoxylin and eosin (H&E) and or IHC, a procedure that does not require particular technological assistance. In the last few decades, the deployment of computer-based, artificial intelligence (AI)-aided imaging analysis to H&E histopathologic slides has increased over time, expanding our capability of image inspection, and improving the identification of image-based tumor classifiers (1–3). Features at both cell and tissue level can be extracted through computational analysis of digitized slides, resulting in a collection of human interpretable features (4, 5). Beyond well-established tumor cell morphologic characteristics, such as nuclear size or number of mitotic cells, novel emerging human interpretable features include spatial arrangement of tissue elements (6), lymphocyte count or spatial distribution (7, 8), and stromal arrangement (9). Promising results are showing the applicability of AI-powered approaches to cancer diagnostics, biomarker development, and prediction of genomic aberrations (1, 8, 9–13).

A different scenario emerges when we consider the development of tools to study the immune and tumor microenvironment (TME). Immune cells critically modulate tumor growth and progression in several ways that can vary according to tumor type, stage of disease, and therapeutic treatment (14–17). Our understanding of the TME has been dramatically improved by multidimensional platforms [e.g., CODEX, imaging mass cytometry (IMC), multiplex IHC; ref. 18], allowing us to detect multiple markers on the same slide and appreciate a previously unexplored immune heterogeneity. One of the key elements provided by such approaches is the localization of immune cells in the spatial context of tumor tissue, and distance metrics reflecting their potential functional interactions. A wide range of novel spatial-related immune features has been emerging recently, including spatial distance, local density measures/clustering and quantitative distribution in tissues, and they may gain even more importance when analyzed in relation to more comprehensive analysis of tumor cell organization.

A complete understanding of the tumor microenvironment requires the integration of H&E analysis and hi-plex approaches. Critical barriers that need to be considered toward an optimal integration of H&E and hi-plex approaches include: (i) information acquired on whole-slide images (WSI) stained with H&E needs to be considered for selection of relevant regions of interest (ROI) required by hi-plex approaches, (ii) H&E is more common in the clinical diagnostic setting, while multiparametric approaches have been so far confined to research studies, and (iii) considerable costs of hi-plex approaches that reduce their applicability in the clinical setting.

Lung cancer is the leading cause of cancer-related deaths in both men and women (19, 20). Despite this bleak scenario, our understanding of the disease is remarkably increasing, screening programs have improved diagnosis, and therapeutic options are encouraging. The most frequent lung cancer type, non–small cell lung cancer (NSCLC), is commonly classified into adenocarcinoma and squamous cell carcinoma based on the cell of origin, but further distinction into histologic patterns has been proposed. For instance, according to histologic examination of architectural features, lung adenocarcinoma is further sorted into lepidic, acinar, papillary, micropapillary, and solid (21), and significant associations of the distinct patterns with prognosis have been documented (22). Of note, more than 80% of the tumors present with a mixture of different histologic patterns (21, 23), making the classification of such patients difficult. In recent years, AI-based approaches have been implemented to sort out the inherent subjectivity of the diagnostic procedure in those cases presenting mixed histologic subtypes (23). Concomitantly, high-resolution analyses aimed at refining the classification of the TME in NSCLC have brought to light a wide range in the extent of immune infiltration, in both lung tumors with different histologies and different regions of the same tumor (24). Efforts are ongoing into the direction of integrating the data collected through multiple platforms into a classification system (25). Recently, a TME score encompassing tumor and immune features has confirmed a significant association between tumor mutational burden and T-cell infiltration in NSCLC and between immune infiltration and survival (26) or response to immunotherapy (8). Of note, tumor features explored by computational imaging analysis of NSCLC specimens have been successfully integrated with transcriptomic information related to the immune infiltration (4).

To improve our capability to describe the NSCLC ecosystem, in this work we aimed at integrating tumor classifiers (generated by AI-powered H&E) and immune classifiers (obtained by multiparametric IMC or AI-powered H&E). By integrating the two approaches, we classified tumor specimens according to tumor cell spatial distribution and analyzed their distinct immune infiltration.

The study population included two cohorts. Cohort 1 included 108 patients diagnosed with NSCLC, whose tumor specimen was spotted on a tissue microarray (TMA; LC2162–US Biomax), containing 209 cores. Patient demographic and histopathologic data are listed in Supplementary Table S1. Of the total spots of the TMA, IMC was performed on the 158 that were localized within the region of ablation. Only spots (n = 116) that passed a staining quality check (including tissue integrity, absence of necrosis, unspecific staining) were considered for downstream analyses. Cohort 2 included 50 whole-tissue slides from 50 patients with NSCLC (one slide/patient) surgically resected at Humanitas Research Hospital between 2015 and 2018. Patient demographic and histopathologic data are listed in Supplementary Table S2. The patients had not received radiotherapy, chemotherapy or immunotherapy before the resection. Written informed consent was obtained from each patient included in the study. The study protocol was in accordance with the ethical guidelines established in the 1975 Declaration of Helsinki and compliant to the procedures of the local ethical committee (protocol no. 63/20).

Four-micron–thick formalin-fixed, paraffin embedded (FFPE) sections from NSCLC tissue blocks or from the TMA were deparaffinized in xylene and rehydrated through a graded alcohol series. The tissue sections were stained with hematoxylin (Histo-Line) for 10 minutes and with eosin (Histo-Line) for 7 minutes. Stained slides were dehydrated through a graded alcohol series before mounting with Eukitt (Sigma-Aldrich). Whole–slide digital scans were acquired using a ZEISS Axio Scan Z1 Slide Scanner at ×20 magnification. Tissues were carefully checked for integrity before further analyses.

The image analysis software QuPath (27) was used to identify tumor cells and to automatically extract their point pattern distribution from digitized H&E images. In cohort 1 (TMA), the entire area of each spot was considered for the analysis. In cohort 2 (i.e., whole slides), for each slide only the tumor region identified by the QuPath's Simple tissue detection tool was considered for the analysis. Cell segmentation was performed with Stardist (28), a public, pretrained deep learning model available as an extension in the QuPath software. To discriminate the tumor cells from the other cells in the tissue, we trained a cell classifier based on the manual annotation of 5–10 stromal or tumor regions containing about 50 to 150 stromal or tumor cells, respectively, adopting the Random trees algorithm and the Train Object Classifier tool embedded in the QuPath software (Supplementary Fig. S1A; ref. 27). The algorithm was fed with more than 50 features per cell [shape features (i.e., area, length, circularity) and intensity features (H&E)]. Accuracy of the classifier in detecting tumor cells was confirmed by staining the consecutive slide with p40 (BIOCARE MEDICAL, catalog no. ACI 3066 A, C, RRID:AB_2858274) and thyroid transcription factor-1 (TTF-1; Ventana Medical Systems, catalog no. 790–4398, RRID:AB_2335981) antibodies (29). To cope with the presence, in our cohorts, of different histologic subtypes (e.g., basaloid, papillary, mucinous, lepidic) and obtain a more accurate image analysis, we generated 10 distinct classifiers according to the tissue and cell morphology. For each image, the classification results were first visually checked, and the accuracy of the AI-predicted tumor cell count was then evaluated using a frame-based comparison method (10). Briefly, the counts obtained by the AI model in 20 randomly chosen patches (100×100 μm) from 10 TMA spots were compared with the manual count from an expert pathologist and the correlation estimated by Spearman correlation analysis (Supplementary Fig. S1B). To quantify the point pattern distribution of tumor cells (expressed in x-y coordinates), for each sample we computed the Ripley's K function, normalized and centered as proposed by Lagache and colleagues (30), taking advantage of the spatstat R package.

The QuPath image analysis software (27) was used to identify T cells from digitized H&E images of Cohort 2. After cell segmentation (described in the previous section), the Points tool embedded in QuPath was used to annotate about 25–30 cells as “T cells” on a single H&E-stained slide, taking as reference its consecutive slide IHC stained with anti-CD3 antibody (Thermo Fisher Scientific, catalog no. MA1–90582, RRID:AB_1956722). The resulting annotations were used to train a random tree–based classifier able to recognize T cells adopting the Train object classifier tool embedded in QuPath. Once trained, the classifier was applied to the whole Cohort 2, and the accuracy of the AI-predicted tumor cell counts was evaluated using a frame-based comparison method (10), where digital counts in 20 randomly chosen patches 200×200 μm from four whole slides were compared with a manual count from an expert pathologist using Spearman correlation (Supplementary Fig. S1C). Heterogeneity across sectioning levels was tested in five different slides from the same tissue block (n = 3; interval 100 μm; Supplementary Fig. S1D). We then computed for each slide (i.e., for each patient) the density of lymphocytes as the number of cells divided by the area of the region within both the tumor and the peritumor, that is, the tissue surrounding the tumor annotation. The combined score was used to quantify the grade of immune infiltration of the tissue and classify the patients’ slides as desert, excluded, or inflamed. Briefly, if the density of lymphocytes within both tumor and peritumor was lower than the respective critical values, the sample was defined as desert; if the density within the tumor was lower than the third quartile, but the density within the peritumor was higher than the median value, the sample was classified as excluded; if the density within the tumor was higher than the third quartile, the sample was inflamed. The region corresponding to the invasive margin was considered as part of the tumor core.

Unconjugated antibodies were conjugated to lanthanide isotypes (Supplementary Table S3) using the MaxPar X8 Antibody Labeling Kits (Standard Biotools) according to the manufacturer's protocol. After conjugation, all antibodies were resuspended in PBS^+/+ with 0.05 mol/L NaN₃ to reach the final concentration of 0.5 mg/mL. For each metal-tagged antibody, the optimal staining condition was determined in a titration assay.

TMA specimens (cohort 1) and 9 WSIs (of cohort 2) were deparaffinized in xylene and rehydrated through a graded alcohol series. Antigen retrieval step was performed in a water bath at 98°C for 20 minutes in Target Retrieval Solution, pH 9 (Agilent), followed by cooling down in distilled water. To block nonspecific binding, the slide was incubated with PBS containing 0.05% Tween, 0.1% Triton X-100, 3% BSA (Sigma), 5% normal mouse serum (Sigma-Aldrich), 5% normal rabbit serum (Sigma-Aldrich), and 5% normal goat serum (Sigma-Aldrich) for 45 minutes. The slide was then incubated overnight at 4°C with a metal-tagged antibodies panel (Supplementary Table S3) and then washed twice in 0.2% Triton X-100 (Sigma) PBS^+/+ for 8 minutes and in PBS^+/+ for 8 minutes. Finally, the TMA was incubated with 0.63 μg/mL Cell-ID Intercalator-Ir (Standard Biotools) for 30 minutes at room temperature washed with distilled water and air dried.

IMC was performed on a Hyperion Imaging System (Standard Biotools). Before laser ablation, a 20× objective lens panoramic optical picture was taken. For the NSCLC tissue microarray, 116 regions of interest (ROI) with an area of 1 mm² were identified to correspond to each single spot of the TMA within the ablation field. For the WSIs, 3 ROIs with an area of 1 mm² were randomly identified within the tumor (consecutive sections stained with H&E were used to distinguish the tumor region from the peritumor area). TMA cores were laser-ablated with a laser frequency of 200 Hz. Data were exported as MCD files and visualized using the MCD viewer software (Standard Biotools).

Count data (MCD files) were converted to TIFF format and segmented into single cells following the IMC analysis pipeline available at https://github.com/BodenmillerGroup/ImcSegmentationPipeline (31). Briefly, the Ilastik v.1.4.0 software was used to interactively train a supervised random forest pixel classifier on multiple image crops, by manually drawing labels for nuclei, membrane or cytoplasm, and background, on the basis of the IMC Cell Segmentation Kit (Standard Biotools) staining. We needed about 10–20 annotations within 5 crops of 500×500 μm to reach an accurate segmentation performance. For each image, the trained classifier was used to create a probability segmentation mask that was then imported in CellProfiler v.4.2.1 software to perform the cell segmentation step with the Propagate algorithm (32), whose accuracy in identifying cells has been validated.

Single-cell segmentation masks were overlaid to the corresponding IMC TIFF image data by ImageJ (RRID:SCR_003070) and the resulting segmented images were then imported on the QuPath v.0.3.2 software. For each spot, the tissue was segmented into stromal and tumor areas by training random trees–based pixel classifiers based on the Pan-cytokeratin staining, while the signal of collagen I was used to annotate fibrotic regions within the tissue. About 10–15 annotations with a mean area of 20,000 μm² were sufficient to train each classifier. TMA spots and sample areas affected by disrupted or necrotic tissue were removed by this and subsequent analyses.

Data on the average ion count per marker (corresponding to the estimation of protein abundance) and on spatial localization were extracted for each cell. Only markers that passed a quality check were retained for the downstream analysis (Supplementary Table S3; Supplementary Fig. S2A). Cell clustering was performed stepwise. First, using Cytomap v.1.4.2, the whole dataset was normalized using the z-score normalization method and, to detect and separate rare cell subpopulations, cells were over-clustered into 100 groups by the NN Self Organizing Map algorithm on the basis of 13 selected cell lineage markers (aSMA, vimentin, CD31, CD45, CD68, CD163, CD3, CD8a, CD4, CD20, CD66a, ki67, FoxP3). Each cluster was then annotated according to their phenotype, and those corresponding to the same cell type were merged, resulting in 11 final cell populations (Supplementary Fig. S2B). To do this, we generated the marker expression heat map hierarchically ordering the clusters with the Pearson correlation distance and cut the dendrogram to identify the main phenotypes of interest. Afterwards, the accuracy of the cell type classification was visually validated on the original images, and some clusters had to be reassigned to the correct phenotype. Marker expression heat maps were made with RStudio (v.2022.07.1+554), using the pheatmap package.

Among the cells annotated as “tumor and other cells” in the first phenotyping step, we classified all cells localized in tumor areas that were neither stromal nor immune cells as “tumor cells” and all the cells localized in the stroma that were neither stromal nor immune cells as “others”. To avoid misclassification of cells along the borders of the tumor nests, the phenotyping step was subjected to adjustments. Briefly, cells that were noted during the first computing as “tumor cells” or “others” in the stroma and located at a distance less than 10 μm from the tumor-nest border were considered part of the tumor area and thus reclassified as “tumor cells” (Supplementary Fig. S2C). Expression levels of pan-cytokeratin confirmed the accuracy of the analysis (Supplementary Fig. S2D). Cells annotated as macrophages were reclustered into 10 groups according to the expression of 11 markers (vimentin, CD14, CD16, CD163, CCR4, CD63, CD68, C1Qa, arginase-1, S100A8, HLA-DR). For those markers (CCR4, CD63, C1Qa, arginase-1) with a low signal-to-noise ratio, the noise was subtracted before performing the downstream analyses. Data integration and phenotype reclustering were performed as previously described (z-score normalization, NN Self Organizing Map algorithm, manual setting of the number of clusters).

Different tissue portions can be differently laser-ablated during the Hyperion workflow. To evaluate the batch effect due to potentially different laser ablation, samples were classified as 1, 2, or 3, according to the level of ablation visually estimated (1, poorly; 2, average; or 3, highly ablated, respectively). We then compared the Uniform Manifold Approximation and Projection (UMAP) of the identified cells, classified according to ablation level, with the UMAP showing the clusters found (Supplementary Fig. S2E). Observed technical variation due to tissue ablation was adjusted by applying the ComBat_seq function (sva R package) directly on the raw count matrix (mean expression of each marker for each cell). The corrected matrix was then imported into Cytomap and the clustering was reperformed.

All statistical analyses were conducted using the software GraphPad Prism 9 (RRID:SCR_002798). Significance was assigned at P < 0.05, unless stated otherwise. Specific tests are indicated in the relevant figure legends. Survival analysis was performed considering disease-free survival (DFS) as outcome measure. All time-to-event endpoints were summarized using the Kaplan–Meier method. Differences in these endpoints between groups were examined using a log-rank test.

All raw data generated in this study is available from the corresponding author upon request.

Visual inspection of tissue slides currently represents the best available method to classify NSCLC specimens into histopathologic subtypes (23). To explore more quantitative histopathologic indicators of tumor cell growth, we implemented an AI-powered digital pathology approach to identify tumor cells on H&E-stained slides (Fig. 1A). On a TMA containing 158 NSCLC tumor cores from 83 patients (Supplementary Table S1), single cells were first detected using Stardist, a public pretrained deep neural network–based model (Fig. 1A, top and middle; ref. 28), and tumor cells were then identified among others through a machine learning–based classifier trained de novo in QuPath (Fig. 1A, bottom). To confirm the accuracy of the AI tool in cell detection, we stained a consecutive section with thyroid transcription factor-1 (TTF-1) and p40, two common tumor markers routinely used to diagnose adenocarcinoma and squamous cell carcinoma (SCC), respectively (Fig. 1B; ref. 29). In both adenocarcinoma (Fig. 1B, top) and squamous carcinoma (Fig. 1B, bottom) subtypes the digital tool accurately identified cells stained positive for the tumor markers. Moreover, the tool did not recognize epithelial cells in normal nontumor lung tissue (Fig. 1C, top), nor normal epithelial cells in peritumor lung tissue (Fig. 1C, bottom). The validity of the tool was confirmed across different histologic patterns, including the mucinous, papillary, and basaloid (Fig. 1D).

The AI-tool allowed us to obtain spatial coordinates for every cell. We then aimed at extracting additional tumor features by quantifying the spatial localization of tumor cells with the K function summary statistic (Fig. 2A; ref. 30). Briefly, the K function describes the spatial distribution of a group of points (corresponding, in our case, to the centroid of each tumor cell) as a function of their distance (Fig. 2B). The resulting K function score (resulting from the distribution of the peak values of the normalized and centered K function curves; ref. 30) expresses how clustered a set of points (or cells) is when compared with a reference set of points with a random distribution (theoretical curve; Fig. 2B and C). Each sample was classified according to the K score, as “uniformly distributed” (Fig. 2A, left, gray circle), “poorly clustered” (Fig. 2A, middle, blue circle) or “highly clustered” (Fig. 2A, right, magenta circle). Another interesting metric related to the tumor cell distribution is the cluster's radius, resulting from the distribution of the radius values corresponding to the peak of the K curves and divided by 1.3 (30), expressing the mean size of the cluster (Fig. 2C). Deployment of the K function calculation to the TMA spots revealed that the cell distribution pattern was highly heterogeneous, with only one adenocarcinoma and one SCC classified as uniformly distributed in the tissue, and the other spots classified as clustered (Fig. 2D, left). The range of K score for these spots was quite large, stretching from 4.80 to 325.81 (Fig. 2D, left), as indication of the heterogeneous distribution of tumor cells in the cohort analyzed. As to the cluster's radius, the population was classified as “small radius” or “large radius” (Fig. 2D, right). The two metrics moderately correlated (n = 116, r = 0.2292, P = 0.0133 by Pearson correlation analysis; Fig. 2E). Notably, among SCCs, there was the highest frequency of highly clustered samples (P < 0.001 by chi square; Fig. 2F).

Integration of AI-powered analysis of H&E slides and multiparametric imaging approaches can provide additional features related to the tumor immune ecosystem. To this aim, we took advantage of IMC and stained one consecutive tissue section of the TMA with metal-labeled antibodies (Fig. 3A; Supplementary Table S3). Following the Hyperion workflow, we could quantitatively detect a panel of 24 markers related to tumor cells, tissue architecture, immune cells, and cell activation (Supplementary Fig. S2A). In addition to cancer cells, based on the expression of lineage markers we could identify nine cell populations. Among the immune cells there were macrophages (Mφ), expressing CD68 and CD163, CD8⁺ T cells, CD4⁺ T cells, T regulatory (Treg) cells, B cells, and proliferating Ki67⁺ immune cells (Fig. 3B; Supplementary Fig. S2B). Stromal populations included αSMA⁺ cells, fibroblasts (expressing lower levels of αSMA and high levels of vimentin) and vessels (Fig. 3B). To investigate how the immune ecosystem changes according to the spatial tissue distribution, we segmented each sample into two regions: tumor-nest (Pan-cytokeratin^pos) and stroma (Pan-cytokeratin^neg; Fig. 3C; Supplementary Fig. S2D) and investigated the lung tumor–immune ecosystem within these two regions. The proportion of all the immune populations analyzed (i.e., Mϕ, CD8⁺ T cells, CD4⁺ T cells, Treg, and B cells) significantly differed within the two tissue regions, being the density higher in the stroma than in the tumor (Fig. 3D, left). Macrophages were the most abundant population both in the stroma and in the tumor-nest, accounting for about 50% of the tissue leukocytes in both the areas (Fig. 3D, right). In the lung, several studies have shown the coexistence of ontogenically distinct Mϕ, both resident (alveolar) and infiltrating (33), as well as Mϕ with different morphologies (34), states of activation and polarization (35). To probe whether infiltration of distinct Mϕ populations correlates with distribution of tumor cells, we performed a reclustering of Mϕ (Fig. 3E), based on the positivity and intensity of 11 phenotypic markers (CD68, CD63, HLA-DR, CD14, CD16, CD163, C1Qa, CCR4, vimentin, S100A8, arginase-1), allowing us to identify 20 distinct clusters (Fig. 3E): Mϕ (c1 and c4 grouped together based on their similarity), expressing CD68 and lacking other activation markers, CD68^hi Mϕ (c12 and c19 grouped together), CD14^hi/CD16^hi Mϕ (c10 and c18 grouped together), CD163^hi Mϕ (c3, c5 and c11 grouped together), Arg1^hi Mϕ (c2), Arg1^hi/S100A8^lo Mϕ (c17), representing the smallest portion of Mϕ, CD63^hi Mϕ (c9 and c20 grouped together), C1Qa^hi Mϕ (c7), Vim⁺ Mϕ (c8 and c13 grouped together), HLA-DR⁺ Mϕ (c6 and c14 grouped together), and S100A8⁺ Mϕ (c15 and c16 grouped together; Fig. 3E). While some clusters occupied definite positions on the UMAP, the Mϕ populations expressing high levels of lineage markers (Mϕ, CD68^hi Mϕ, CD14^hi/CD16^hi Mϕ and CD163^hi Mϕ) were closer and overlapping, suggesting that they might represent the same population at different maturation states (Fig. 3F), confirming previous evidence obtained by trajectory analyses of transcriptomic data (36, 37). By UMAP visualization of the cluster arrangement with respect to their tissue location (Fig. 3G), the clusters expressing S100A8 and arginase-1 (c15, c16, c2, c17) populated the tumor-nest and the stroma in a comparable way (Fig. 3G and H), in contrast to the other clusters that were significantly more abundant in the stroma (Supplementary Fig. S2F). This result confirms that macrophages with different profiles occupy different topological positions.

The preferential location of immune cells in the nontumor, stromal region could be due to several factors and integration of the AI-based tumor classifiers and stromal/immune classifiers obtained by IMC could provide additional insights. To assess whether the extracellular matrix impacts on the distribution of immune cells in NSCLC, we refined the tissue annotation by using the collagen I signal, the predominant component of the fibrotic tumor stroma. Collagen distribution was heterogeneous across samples, with some spots almost totally covered by collagen (Fig. 4A, left) and others with very little amount of it (Fig. 4A, right). Classification of samples into highly fibrotic or poorly fibrotic evidenced a higher frequency of immune cells in highly fibrotic samples, in contrast to tumor cells that were more frequent in poorly fibrotic samples (Fig. 4B). The extension of the collagen area was neither correlated to the tumor K score (Fig. 4C, left) nor to the cluster's radius (Fig. 4C, right), suggesting that the way how tumor cells cluster was not related to the presence of a fibrotic stroma. Of note, the percentage of immune cells was higher in highly clustered compared with poorly clustered samples (P = 0.0488 by Mann–Whitney test; Fig. 4D). This difference was primarily due to adaptive cells. In fact, a significantly higher density of CD8⁺ T cells, B cells, and T regulatory cells was found in the stromal regions of highly clustered samples (P < 0.05 by Mann–Whitney test; Fig. 4E and F). Collectively, paired analysis of the tumor and immune ecosystem by AI-powered H&E and IMC suggests that the spatial distribution of cancer cells is significantly related to the ability of adaptive immune cells to infiltrate the tumor.

Histopathologic analysis of WSIs represents the common setting of the clinicodiagnostic routine. IMC deployed to whole slides confirmed the preferential localization of immune populations in the stroma rather than in the tumor regions (Supplementary Fig. S3A). To explore the value of the spatial classifiers identified in the TMA in a format suitable for translational research, we deployed the AI tool to H&E-stained WSIs from a cohort of 50 patients with NSCLC (Supplementary Table S2; Fig. 5A). We then calculated the K function for the WSIs (Fig. 5B and C), which allowed us to classify the samples according to K score and cluster's radius (Fig. 5C, bottom), confirming that translation to the WSI setting is applicable. Because paired analysis of the tumor and immune ecosystem by AI-powered H&E and IMC had pointed to a relevant distribution of T cells according to the tumor distribution (Fig. 4E and F), we aimed at extracting features related to the T-cell compartment exploiting a machine learning classifier to identify T cells in H&E slides. The digital tool was trained to recognize T cells on H&E images by labeling cells based on a set of images obtained from tissues stained with an anti-CD3 antibody (Fig. 5D, left) and further implemented with segmentation of the tissue in tumor and peritumor compartments (Fig. 5D, middle, right). Analytic precision of the tool was evaluated on distant sections within individual tumor blocks (Supplementary Fig. S1D). Collectively, the digital tool predicted the density of T cells on H&E slides and, according to their preferential location in the tumor or peritumor areas, it allowed to classify each specimen as immune “desert”, “excluded” or “inflamed” (Fig. 5E and F), a common immune-based classification that has prognostic value in different tumor types (38). Patients classified as immune inflamed presented with a significantly better prognosis compared with excluded patients (Fig. 5G; P = 0.064 by Mantel–Cox test). This result was confirmed on a larger cohort (Supplementary Fig. S3B). Combination of the tumor classifier with the immune one resulted in a significant better prognosis of patients classified as inflamed/K^hi (i.e., concomitantly having an inflamed and highly clustered tumor) compared with the excluded/K^lo (i.e., concomitantly having an immune excluded and poorly clustered tumor; Fig. 5H; P = 0.046 by Mantel–Cox test). Moreover, in patients classified as inflamed/K^hi there was a significantly lower frequency of recurrence (Fig. 5I). These results confirm that paired H&E and hi-plex approaches can provide composite classifiers with clinical relevance.

Digital pathology and AI-guided approaches have received great attention lately, both in academic medical centers and clinical pathology units, for ameliorating the current workflows and improving on manual histopathologic scoring of cancer tissues (1–3, 39). AI approaches, especially when weakly supervised, provide the opportunity to disclose histologic information that goes beyond visual perception, ultimately outdoing what pathologists do (39). Efforts are ongoing into the direction of integrating AI-powered H&E with hi-plex platforms, most used to analyze the immune ecosystem (40, 41). For instance, hi-plex platforms generally require the selection of regions of interest (ROI) and are more suitable to the analysis of samples spotted on TMAs, while WSIs are the ideal setting for diagnostic-clinical image inspection. In this work we combined standard H&E and hi-plex IMC and by AI-guided tools, we extracted, from NSCLC slides, quantitative and standardized subvisual features that could have been underestimated by eye. Our integrative approach was preliminarily deployed to a TMA and further validated on WSIs.

The introduction of IHC staining detecting the markers TTF-1 and p40 to discriminate distinct tumor subtypes has significantly improved the diagnostic accuracy for NSCLC. A few issues remain though, such as the interpretation of results, especially in tumors presenting with mixed subtypes (29) or loss of the marker in poorly differentiated tumors (42). We exploited AI-powered H&E to detect tumor cells and extract their spatial coordinates, which was then used to extract a tumor topological classifier related to how the cells distribute in the environment. Not only did the digital tool accurately identify tumor cells, it could also improve the analytic output of the diagnostic procedure. Therefore, the leverage of an AI algorithm generating a tumor H&E classifier mitigated some of the challenges associated with manual assessment and improved the accuracy of tumor-cell recognition. Recently, a deep convolutional neural network has been trained to classify WSIs from the Cancer Genome Atlas in the prevalent subtypes of lung cancer (43). Collectively the recognition of tumor cells by algorithms deployed to H&E regardless of the marker that tumor cells express, may aid in those circumstances in which a reliable tumor cell marker is not available. Besides recognizing tumor cells, the K function statistic (44) captured the modality of tumor cell distribution, providing a spatial classification, indicating the degree of tumor cell clustering. Of note, tumor cells in squamous tumors were more clustered, which is in line with the histopathologic features of this NSCLC subtype.

Multidimensional IMC allowed us to capture additional information related to the TME. First feature we caught was the regionality of immune cell localization, namely their preferential localization in the stroma. Immune exclusion, that is, the confinement of immune cells outside the tumor bed (45), is a relevant feature of the TME and is commonly used to classify tumors into cold/hot or inflamed/excluded. Deep-level analysis of immune populations evidenced additional features. A focus on macrophages, by using common myeloid markers (e.g., CD68, CD14, CD16) and some more specific ones (HLA-DR, S100A8, Arg1, C1Qa), allowed us identifying more than 10 populations. Although a parallel with key Mϕ subtypes described in recent transcriptomics analyses could be imperfect, due to inherent differences between gene profiles and protein expression, we confirmed some common elements. Immature (e.g., S100A8⁺ and CD14⁺CD16⁺Mϕ) and mature (CD163⁺ CD68⁺ Mϕ) macrophages coexist in the same tissue (37, 46), and can be diversified on the basis of their localization, e.g., only S100A8⁺ macrophages infiltrated the tumor region. The regional specialization of macrophage subtypes in tumor tissues is an important feature that can control their functional profile and correlate with survival benefit, as shown in melanoma (47) and in human colo-rectal liver metastasis (37). In line with previous studies, a considerable fraction of the Mϕ granularity was accounted by very similar populations (occupying the same position on a UMAP representation), possibly representing the same cell in different maturation states (36, 37). Another interesting aspect emerging from the integrative analysis was that tumor cells were significantly more frequent in poorly fibrotic tumors, while the presence of fibrosis did not correlate with the K score (i.e., the distribution pattern of tumor cells). These results suggest that a fibrotic stroma impacts on tumor cell number but not on their clustering, while it is not the case for immune cells, particularly the adaptive ones, whose density is affected by fibrosis, as shown already in the literature (48).

The next step, that is, the combination of H&E-based tumor classification and IMC-based immune classification resulted into a composite classifier that was tested as prognostic. Patients classified as immune excluded presented with a significantly worse prognosis compared with inflamed patients. The intermediate profile of desert patients could be due to other immune factors not accounted, such as the presence of myeloid cells or other immunosuppressive populations. Prediction models have been increasingly used to infer the number of immune cells in H&E slides (8), which are widely available for virtually every patient and more suitable to weakly supervised deep-learning algorithms. Their performance is not the same for all the cells and they tend to work better for lymphocytes than myeloid cells. We recently explored the efficacy of deep learning approaches in recognizing morphologic features of macrophages in human colorectal liver metastasis sections (13) and we experienced some issues due to the granularity of the staining, separating adjacent cells, and avoid oversegmentation of large cells.

Our integrative workflow provides an example of how the deployment of AI algorithms to common H&E slides and integration with hi-plex approaches can readily transform variables into insights. Although IMC was the approach we used, the workflow described here provides a general framework for the integration of AI-powered H&E with any hi-plex assay. Large-scale validation and correlation with clinical variables may ultimately translate these features into novel classifiers, to predict prognosis, genotype, and response to therapy from images. To date, AI-powered H&E has received FDA approval for diagnosis of prostate cancer (49). In gastrointestinal cancers, microsatellite instability can be predicted directly from H&E-stained histology slides (50, 51), an attractive application, if we consider that the microsatellite instability status determines patient response to immunotherapy (52).

Even though AI-powered H&E seems easily deployable in a clinical setting, its use in the oncology biomarker scoring has been negligible (39, 53), primarily due to the lack of large-scale validation studies encompassing multiple cohorts and tumor types. Our results, once validated in a larger cohort, may improve the design of diagnostic algorithms, performing more accurate patient classification into clinically relevant groups, leveraging spatial biology for precision oncology.

A. Rigamonti reports grants from AIRC Foundation for Cancer Research in Italy during the conduct of the study. A. Santoro reports other support from BMS, Servier, Gilead, Pfizer, Eisai, Bayer, Merck Sharp & Dohme, Sanofi, Incyte, and Takeda, Roche, AbbVie, Amgen, Celgene, AstraZeneca, Arqule, Lilly, Sandozd, and Novartis outside the submitted work. A. Mantovani reports personal fees from Ventana, Pierre Fabre, Verily, Abbvie, AstraZeneca, Verseau Therapeutics, Myeloid Therapeutics, Third Rock Venture, Imcheck Therapeutics, Ellipses, Novartis, Roche, Macrophage Pharma, Biovelocita, Merck, Principia, BioLegend, Olatec Therapeutics, Moderna, Henlius; grants from Novartis; other support from Cedarlane, Hycult, eBioscience, BioLegend, Abcam, Novus Biologicals, Enzo Life, and Affymetrix outside the submitted work; in addition, A. Mantovani has a patent for WO2019057780 “Anti-human migration stimulating factor (MSF) and uses thereof” pending and issued, a patent for WO2019081591 “NK or T cells and uses thereof” pending and issued, a patent for WO2020127471 “Use of SAP for the treatment of Euromycetes fungi infections” pending and issued, and a patent for EP20182181.6 “PTX3 as a prognostic marker in Covid-19” pending and licensed to Diasorin. F. Marchesi reports grants from Associazione Italiana per la ricerca sul cancro (AIRC) and Italian Ministry of Health during the conduct of the study. No disclosures were reported by the other authors.

A. Rigamonti: Formal analysis, investigation, methodology, writing–original draft. M. Viatore: Formal analysis, validation, methodology, writing–original draft, writing–review and editing. R. Polidori: Data curation, formal analysis, validation, methodology, writing–original draft, writing–review and editing. D. Rahal: Resources, methodology, writing–review and editing. M. Erreni: Methodology, writing–review and editing. M.R. Fumagalli: Methodology. D. Zanini: Methodology. A. Doni: Methodology. A.R. Putignano: Methodology, writing–review and editing. P. Bossi: Resources, data curation. E. Voulaz: Resources. M. Alloisio: Resources. S. Rossi: Resources. P.A. Zucali: Resources. A. Santoro: Resources. V. Balzano: Resources. P. Nisticò: Resources. F. Feuerhake: Supervision, writing–review and editing. A. Mantovani: Supervision. M. Locati: Supervision, writing–review and editing. F. Marchesi: Conceptualization, supervision, funding acquisition, writing–original draft, project administration, writing–review and editing.

This work was supported by Fondi Ministero della Salute – Progetti di Rete ACC RCR-2020–23670066 (to A. Mantovani and P. Nisticò); Associazione Italiana per la ricerca sul cancro (AIRC), under grant AIRC IG 2020 ID. 24393 (to F. Marchesi), under fellowship ID. 25290 2020 (to A. Rigamonti), and under AIRC 5×1000 21147 (to A. Mantovani); the Italian Ministry of Health: RF-2019–12369142 (to F. Marchesi). The funding agency had no role in design of the study, collection, and analysis of data. We thank Dr. Fabio Giavazzi for help in computational analyses, and Dr. Fabio Pasqualini and Fabio Grizzi for assistance with the AxioScan equipment. The PhD student R. Polidori was supported by the PhD program in Experimental Medicine of the University of Milan.