Origins of image analysis and multiplexed imaging techniques

Tissue biomarkers have been employed in biomedical research and in the clinic for a long time. Biomarkers can be broadly classified as prognostic or predictive. Prognostic biomarkers provide relevant tumor behavior information about the probability of survival regardless of the treatment received. On the other hand, predictive biomarkers are able to predict outcomes in a treatment-specific fashion and are therefore very helpful for deciding patient management. Tumor morphology and histological classifications are indeed the most basic tissue biomarkers. As immunohistochemical and immunofluorescence (IF) methods have advanced, it has become possible to discern cell populations and measure their densities within tissue sections. Immunohistochemistry (IHC) has allowed the quantification of immune cell populations and immunomodulatory molecules that have prognostic and predictive power. It is, therefore, important to measure biomarkers in an analytically valid and reproducible manner.

For decades, semiquantitative scoring systems based on pathologists’ estimations using light microscopy have served as the conventional approach for IHC analysis. Moreover, these evaluator-dependent scoring methods are currently the gold standard for assessing several clinically relevant biomarkers, such as PD-L1 expression.11 In addition to other preanalytical and analytical factors,12 it is known that the “manual scoring” by pathologists is subjective, has limited reproducibility and lacks a true quantitative output.13 Therefore, there have been efforts to develop analytical quantification methods.14 15 Although several of these methods can evaluate relevant clinical biomarkers,16–18 they have not been widely adopted in clinical practice.

No single biomarker has demonstrated absolute prognostic or predictive power, which supports the need to integrate data from multiple markers and/or nonredundant metrics. Therefore, simultaneous measurements of diverse markers in a single slide and spatial associations among cell subpopulations could improve biomarker discovery. In addition to the single marker IHC and IF, novel techniques enable researchers to investigate the spatial architecture of the TIME in a more comprehensive manner. Multiplexed and hyperplexed IHC-based,19 20 IF-based,21–23 mass spectrometry-based24 and spatial transcriptomics/proteomics-based25–27 technologies increase the number of simultaneously identified proteins or mRNA transcripts, and therefore, the number of potentially quantified immune and cancer cell populations (figure 1). In addition, the incorporation of functional markers allows the interrogation of specific cell states or biological processes within the cell populations of interest, increasing the significance and dimensionality of the data.

• Download figure
• Open in new tab
• Download powerpoint

Figure 1

Top: Overview of the design, performance and analysis of a tissue imaging experiment to analyze the tissue immune microenvironment (TIME). Bottom: Table summarizing the main technologies used to perform spatial TIME analyses. Figure partially created with BioRender.com. The costs are approximately defined as follows: low<US$500/sample; mid=US$500–US$1500/sample; and high>US$1500/sample. FFPE, formalin-fixed paraffin-embedded; FFT, fresh frozen tissue; IF, immunofluorescence; QC, quality control.

However, these technologies are associated with challenges related to immune marker selection, assay validation, protocol optimization and analytical platform selection. First, multiplexed assays require complex development, control and validation processes prior to implementation.28 Preanalytical considerations, including time from tissue collection to tissue fixation (cold ischemia time), fixation and tissue processing methods, should be carefully planned. In addition, the storage conditions and age of FFPE tissue blocks can influence the antigenicity of target markers.29 30 Second, the scientific question to be addressed should guide the multiplexed platform and the panels of markers to be analyzed. Platforms differ in their relevant technical characteristics, such as spatial resolution, multiplexing capacity, and workaround time. Furthermore, the source of tissue samples dramatically influences the results. For example, a tissue microarray allows the analysis of different samples in a single FFPE tissue section, ameliorating technical variability issues. On the other hand, whole slide images can include larger tissue areas and enable a deeper understanding of intratumor spatial heterogeneity. Last, different analysis pipelines will determine the computing power requirements. A standard multiplex marker assessment pipeline is shown in figure 1.

Single-plex and multiplex marker quantification analysis pipeline

Tissue/cell segmentation and quantification of immune populations are generally the initial steps in most multiplex tissue-based analysis pipelines. Different analysis strategies have unique advantages and are subject to different limitations. Marker quantification relative to the measured area (marker expression per square millimeter), cell density (cell per tissue area such as square millimeter) or the fraction of cells relative to the total number of analyzed cells (proportion of cells) are commonly used metrics. The quantification of individual immune cell populations can be used to calculate their prognostic and predictive potential, to identify biological differences among tumor types or to analyze their variation longitudinally over time in the same tumor/patient after therapeutic interventions. In addition, information on the global infiltration of various immune populations is needed to design complex spatial analysis strategies.

Pixel-based marker quantification

Digital images are numeric arrays of pixels that contain information on the cells present in the tissue. The pixel information present in these arrays will vary according to variations in cell type composition. Therefore, quantifying normalized pixel intensities without considering individual cell units is a commonly employed and effective strategy for analyzing multiplex tissue images. Pixel-based analyses usually begin with image preprocessing steps that are followed by marker normalization to ensure intraslide and interslide comparability. Next, different tissue compartments or areas are defined according to the pixel intensity of specific markers. For example, cytokeratin can be used to define the pixels belonging to the cancer cell compartment, vimentin can be used to identify the surrounding stromal area, or CD3 can be used to delineate the T-cell compartment. Compartment definition usually requires setting adequate marker thresholds and may also require one or several rounds of image preprocessing or filtering to accommodate specific compartment shapes, smooth out compartment edges and remove any artifactual holes or defects in the compartment mask. Subsequently, the pixel intensity of the marker of interest is averaged across the selected tissue compartment and divided by the total compartment area. This metric captures both the number of positive cells/structures that express the marker and their intensity in positive areas. Complex colocalization hierarchies can be used to analyze specific colocalized markers. The main advantages of pixel-based quantification are its simplicity, high precision, objective nature, and reproducibility, in addition to its simple computational requirements. In addition, pixel-based scoring is not restricted to intracellular markers and is well suited for quantifying extracellular immune mediators, such as cytokines and other relevant macromolecules, and complex biological structures, such as neutrophil extracellular traps (NETs), which are webs of extracellular nuclear DNA associated with proteins extruded by dying neutrophils.31 However, pixel-based approaches that quantify immune markers using average fluorescence intensity scores that are not standardized, lack the interpretability of cell-based quantification and are affected by the size of the tissue compartment being analyzed. These drawbacks limit its use in subsequent spatial cell‒cell interaction analyses (see table in figure 2). By using pixel-based quantitative IF (QIF) platforms, such as automated quantitative analysis (AQUA),32 we have investigated the immune microenvironment of a variety of tumors.33 34 For example, pixel-based QIF was used to analyze various immune populations in glioblastomas treated with immunotherapy. Changes in the mean fluorescence intensity of several immune markers were noted.33

• Download figure
• Open in new tab
• Download powerpoint

Figure 2

Overview of the two main strategies (pixel-based and cell-based) used to quantify immunofluorescence markers. Top: Flow chart comparing the analytical pipelines. Bottom: Table displaying the feature comparison between methods.

Cell-based marker quantification

Cells are considered the basic biological unit of tissue and generally consist of three identifiable structural components: the cell membrane, the cytoplasm, and the nucleus. Image analysis algorithms can be applied to segment the membrane and nuclear compartments of cells, which allows single-cell high-dimensional morphology analysis. In contrast to pixel-based quantification, after image preprocessing steps, cell segmentation processes are carried out before proceeding with the analysis. Although still in use, algorithms based on watershed segmentation,35 36 active contour segmentation,37 graph cut segmentation38 or Voronoi tessellation39 (among others) have recently been superseded by deep learning-based segmentation models.40 Most of these algorithms rely on nuclear markers such as DAPI (4',6-diamidino-2-phenylindole) or Hoechst stain to perform nuclear segmentation. Each detected nucleus is then expanded to provide an approximation of the full cell area, constrained by the distance and proximity of neighboring cells. A list of several features derived from cell morphometry and intensity measurements is then calculated for each cell. Alternatively, segmentation methods that use nuclear and membrane markers to feed deep learning-based classifiers, such as the DeepCell algorithm, have been developed as powerful cell segmentation tools.41 They run pretrained networks (eg, Mesmer) to generate cell masks.41 Quality checks of the cell outlines can be visualized by overlaying the cell masks over the original images.

After cell segmentation and feature extraction, the cells are labeled according to their marker expression profile. Expression-matrix approaches rely on the annotation of the average intensity of different markers for each of the cells in the sample. Intensity thresholds can then be specified manually by researchers or by distribution histogram-based methods, such as the Otsu method or Renyi entropy.42 43 Other alternatives not requiring thresholding include various clustering techniques (hierarchical, density-based, and graph-based techniques, such as the Louvain or Leiden algorithms). Alternatively, user-supervised, image-based machine learning algorithms such as CellSighter44 can be trained to classify cells in experiments. These algorithms work directly with images, considering other morphometric variables in addition to marker intensity in the cell segmentation process. After major cell lineage identification, proteins of interest not used for the phenotyping process can be quantitatively scored within the defined cell population (such as granzyme B, PD-1 or PD-L1).

The main advantage of cell-based approaches is that they provide biologically interpretable metrics. In addition, these methodologies allow subsequent complex cell‒cell interaction analysis pipelines. Cell-based pipelines are the mainstays of commonly employed multiplex IF analysis software, such as QuPath,45 Halo (Indica Labs), Phenoplex (Visiopharm) and InForm (Akoya Biosciences). However, the evaluation of individual cell populations using tissue segmentation strategies is subject to prominent associated challenges. These include difficulties in systematically obtaining accurate single-cell identification/phenotyping across tissues and samples, the use of extremely labor-intensive morphology-assisted quality control steps to ascertain adequate performance, the presence of bias secondary to dissimilar marker expression levels/patterns and cell size or morphology, the loss of valuable marker information in the thresholding process through binarization, and the oversight of information present between individual cells in interstitial/extracellular areas that may have biological relevance. Moreover, cell segmentation strategies have associated computational costs. As a result of increased complexity and variability in results between methods, reproducibility may be an issue.

Cell-based quantification has grown in popularity in recent years. For example, a cell-based strategy has been used to quantify CD8+T cells in various tumor types and to correlate cell density with the presence of NETs.46

Quantifying spatial intratumor heterogeneity

Tumors are subjected to evolutionary pressure that leads to intertumor and intratumor heterogeneity (ITH).47 In recent years, ITH has received increasing attention due to its association with cancer progression and response to treatment, including its response to checkpoint inhibitors.48–52 ITH quantification is a field of active research. From a genetic perspective, multiple methods have been developed to quantify different tumor subclones as a proxy of ITH.53–56 In addition, efforts have been made to analyze ITH at the transcriptomic level.57 58 At the spatial level, tumors also display heterogeneity. However, analyzing spatial-ITH (s-ITH) is a complex task. At the spatial level, tumors can deviate from perfect homogeneity in several different ways (figure 3), implying that measuring heterogeneity requires the integration of several different analytical strategies. Furthermore, s-ITH refers not only to cancer cells but also to multiple types of immune and nonimmune cells (such as fibroblasts and endothelial cells). In addition to variations in intrinsic cell properties and cell density, s-ITH should also reflect the diversity of spatial configurations present in the sample.

• Download figure
• Open in new tab
• Download powerpoint

Figure 3

Spatial intratumor heterogeneity analysis (s-ITH). Top: (A) Spatial heterogeneity can display different patterns. Some patterns are associated with high global heterogeneity while others are not. (B) Flow chart with the main outline of an s-ITH analytical pipeline. Bottom: Table describing different commonly used s-ITH methods. Figure partially created with BioRender.com.

Broadly, s-ITH can be analyzed globally (considering a tumor as a whole population) or locally (considering local spatial communities within the tumor). Although cells are the basic unit of ITH analysis at the genomic or transcriptomic level, no consensus exists on the spatial unit of measure used to perform local s-ITH. Most local s-ITH approaches are based on the selection of different areas by the researcher or the segmentation of the image into areas of constant size (image tiling or tessellation). Relevant metrics that capture the density or proportion of various cellular phenotypes are estimated in each of these areas. Metrics to measure s-ITH differ in their complexity and robustness. The metric choice can heavily influence the data structure, properties, and results of the analysis. ITH metrics can be further divided into statistical, image-based, and graph-based metrics. Statistical methods represent some of the simplest ITH metrics, including dispersion estimators such as variance, standard deviation (SD), range, and coefficient of variation. More intricate methods, such as entropy metrics59 or mixed effects models,60 can also be applied. Image-based methods analyze the spatial distribution of single markers to approach the level of spatial homogeneity or clustering. Textural image analysis61 and spatial autocorrelation metrics62 are examples of image-based approaches (see table in figure 3). Other approaches represent tissue images as graphs, where cells are nodes and edges are cell distances. Each cell or node is assigned to a specific cell type. Graph modularity scores measure how different cell types are spatially arranged into homogeneous or heterogeneous communities.63

Excluding entropy and graph-based metrics, all other s-ITH metrics are restricted to analyzing single cell types. In addition, heterogeneity measurement is heavily influenced by the analyzed tissue area and the size and number of tiles used to implement the analysis. Furthermore, the number of cell phenotypes and their density also impact the metric results.

These types of analyses revealed the impact of s-ITH on survival in patients treated with immunotherapy. For example, Lopez de Rodas et al combined a local tessellation approach to calculate a global per-sample Rao quadratic index (an entropy metric) to classify non-small cell lung cancer (NSCLC) samples into high-heterogeneity and low-heterogeneity groups. The authors demonstrated that patients with tumors displaying high heterogeneity had worse survival after treatment with checkpoint inhibitors.64

Simple spatial cell‒cell interaction

Compared with other multiplexed technologies, image-based multiplexed platforms allow direct interrogation of cell‒cell spatial mapping. Cell–cell spatial interactions are required for many immunologically relevant processes, such as cell migration, antigen presentation and cell cytotoxicity. Spatial Cartesian mapping can reveal biologically and clinically relevant complex immune architectures that cannot be identified by cell population quantification.

Most spatial interaction analysis pipelines start by measuring the distance between all the cells from two populations of interest, obtaining a cell-to-cell distance matrix. Basic spatial interactions include calculating the average distance, the average minimum distance, or the average maximum distance between two cell populations (figure 4). Other approaches require the user to define biologically relevant interaction distances and select a reference and target cell. Afterward, the number of target cells within the defined distance from every cell of origin is calculated. Alternatively, other approaches, such as the G-cross function,65 infer the probability of encountering a target cell across the distance radius. Notably, most Euclidean cell‒cell distance metrics obtained from segmented tissue samples are highly dependent (and inversely correlated) on the density of the cells being analyzed. Thus, careful analysis using correlation functions or multivariate models should be incorporated to ascertain their contribution beyond cell counts. Alternatively, to account for the influence of the total cell density on spatial correlation metrics, several methods have been developed. Ripley’s K function66 (or similar variants including L, F and G functions) can be used to calculate the colocalization of one or several cell populations within the distance radius K compared with a complete spatial randomness distribution. This random background is usually obtained from a Poisson process, by homogeneous distribution sampling or by iteratively permuting cell labels. In addition, random distributions and colocalization metrics can be calculated considering tissue hierarchies (eg, tumor and stromal compartments).67 One of the main limitations of these approaches is that the colocalization distance threshold must be selected. This is usually an arbitrary choice, as only rarely is there enough evidence to support an explicit biologically relevant distance threshold. To address this limitation, numerous radius sizes can be interrogated using cross-K functions, as implemented in the R package SPIAT.68 Moreover, distance analysis is computationally intensive and is usually performed for restricted cell populations. Despite being a powerful technique, analyses restricted to two populations may not be sufficient to capture the complexity of spatial interaction architectures.

• Download figure
• Open in new tab
• Download powerpoint

Figure 4

Simple cell spatial interaction analysis and complex (neighborhood) interaction analysis. Top left: Table showing the different methods used to analyze cell-to-cell spatial interactions. Top right: Explanatory diagrams of simple cell-to-cell spatial interaction analysis. Bottom: Explanatory diagram of closest neighbors and UTAG neighborhood analysis approaches. UTAG, unsupervised discovery of tissue architecture with graphs.

For example, using radius interaction analysis, it has been shown that close contacts and interactions (defined as having a radius of less than 45 µm) between PD-1+CD8+T cells and PD-L1+cancer cells are associated with a poor prognosis in patients with endometrial carcinomas.69

Multiple spatial interactions: cellular neighborhoods

The variety of cell types present in the TIME requires the visualization and quantification of several immune and cancer cell markers in a single tissue section. Many immune-relevant functions within the TIME are organized into highly complex tridimensional structures that are characterized by a conserved composition and spatial distribution of immune cells. These conserved, distinct structures can be described as two-dimensional (2D) representations using cellular neighborhoods. For example, the organization of several B and T-cell lymphocyte populations as well as dendritic cells into tertiary lymphoid structures is associated with favorable prognosis in several cancer types.70 71 Systematic spatial mapping of various cell types in a nonrandom fashion is the basis of cellular neighborhood analysis.

Broadly, cellular neighborhood analysis strategies typically incorporate machine or deep learning strategies and can be classified as unsupervised or weakly supervised. Unsupervised approaches do not take into account image metadata to identify tissue neighborhoods. Among these algorithms, one of the most popular unsupervised neighborhood profiling algorithms requires the annotation of the closest neighbors surrounding every single cell in the sample. These closest neighbor patterns are then used to feed a clustering algorithm to find reproducible spatial colocalization patterns.72 This approach has been used to decompose the immune microenvironment in colorectal carcinomas using the codetection by indexing (CODEX)73 labeling technique as well as to characterize complex cellular neighborhoods in NSCLC.74 Other algorithms, such as unsupervised discovery of tissue architecture with graphs (UTAG),75 are based on adjacency and feature matrices. In this approach, a distance threshold to consider two cells as spatially related must be defined by the researcher. UTAG allows the discovery of ‘tissue microanatomical domains’ that can be used to decompose tissue into meaningful spatial interaction domains. Other approaches inspired by natural language processing, such as SpatialLDA,76 analyze the spatial neighbors of cells to infer different spatial microenvironments (called “topics”), just as groups of words are usually associated with specific topics in a text document. This approach has been applied to identify novel spatially recurrent cell‒cell microenvironments associated with recurrence in colorectal cancer.77

Approaches such as that implemented in CytoMAP78 analyze tumor neighborhoods in a grid-like manner. The cell composition and marker intensity were calculated for all the tiles in the sample, and the results were clustered using self-organizing maps.

Recently, weakly supervised deep learning methods, such as NaroNet,79 have been implemented to perform cellular neighborhood analysis. This trend is driven by image metadata (including clinical information, molecular data, and histologic data) used as training labels to obtain biologically relevant spatial patterns. These patterns are specific for the selected labels. In contrast to other methods that consider a single definition of neighborhood, NaroNet’s spatial interaction networks are arranged in hierarchies of increasing size and complexity (local phenotypes, neighborhoods, and areas). This approach allows a broader analysis of complex cell‒cell interactions. Using NaroNet, complex cell‒cell spatial interaction patterns associated with endometrial cancer relapse have been analyzed,80 identifying that CD8+T cells, Tregs and macrophage-rich neighborhoods are associated with better clinical outcomes.

The choice of the neighborhood approach usually relies on the scientific hypothesis to be addressed. Unsupervised approaches are good options for pattern discovery studies while weakly supervised strategies can be employed to increase the clinical relevance of results. The notion that neighborhood analysis is complex and computationally intensive should also be considered. In addition, all algorithms require the adjustment of hyperparameters, which may be performed in a user predefined manner or through cross-validation methods. These parameters influence the ability of the algorithm to identify different spatial patterns present in the tissue. Although less prone to bias, hyperparameter cross-validation increases the computational costs of the pipeline. However, flawed hyperparameter adjustment may lead to reproducibility issues. A major limitation of computationally complex neighborhood analysis relies on the numerous data required to train stable and reproducible models. Neighborhood identification methods that rely on deep learning may require an exceedingly large number of image examples to be trained and validated. Although this issue may partly be alleviated with the use of synthetic images created through several different methods,68 81 it still represents a major limitation to its widespread implementation. In addition, there is a lack of uniformity in the computational rules and hyperparameters necessary to capture biological or clinically relevant information. These models are also generally subject to the same limitations of single-cell segmentation analysis described above. To overcome some of these limitations, a recent model called ‘Graph of Graphs’ has incorporated a manifold learning approach to achieve accurate classification of tumors from patients with lung cancer sensitive and resistant to PD-1 axis blockers. The authors used a segmentation-independent strategy in which tissue microdomains or “patches” from tumor slides stained with a 29-marker imaging mass cytometry panel were used to generate multimarker correlation vectors that represent the local tissue organization.82