SPARTIN: a Bayesian method for the quantification and characterization of cell type interactions in spatial pathology data

Introduction: The acquisition of high-resolution digital pathology imaging data has sparked the development of methods to extract context-specific features from such complex data. In the context of cancer, this has led to increased exploration of the tumor microenvironment with respect to the presence and spatial composition of immune cells. Spatial statistical modeling of the immune microenvironment may yield insights into the role played by the immune system in the natural development of cancer as well as downstream therapeutic interventions.

Methods: In this paper, we present SPatial Analysis of paRtitioned Tumor-Immune imagiNg (SPARTIN), a Bayesian method for the spatial quantification of immune cell infiltration from pathology images. SPARTIN uses Bayesian point processes to characterize a novel measure of local tumor-immune cell interaction, Cell Type Interaction Probability (CTIP). CTIP allows rigorous incorporation of uncertainty and is highly interpretable, both within and across biopsies, and can be used to assess associations with genomic and clinical features.

Results: Through simulations, we show SPARTIN can accurately distinguish various patterns of cellular interactions as compared to existing methods. Using SPARTIN, we characterized the local spatial immune cell infiltration within and across 335 melanoma biopsies and evaluated their association with genomic, phenotypic, and clinical outcomes. We found that CTIP was significantly (negatively) associated with deconvolved immune cell prevalence scores including CD8+ T-Cells and Natural Killer cells. Furthermore, average CTIP scores differed significantly across previously established transcriptomic classes and significantly associated with survival outcomes.

Discussion: SPARTIN provides a general framework for investigating spatial cellular interactions in high-resolution digital histopathology imaging data and its associations with patient level characteristics. The results of our analysis have potential implications relevant to both treatment and prognosis in the context of Skin Cutaneous Melanoma. The R-package for SPARTIN is available at https://github.com/bayesrx/SPARTIN along with a visualization tool for the images and results at: https://nateosher.github.io/SPARTIN.

1 Introduction

While the staining and examination of tissue samples has been a ubiquitous medical practice for decades, it has only been relatively recently that high-resolution pathology images of such stainings have begun to be explored through the lens of formal quantitative models (Heindl et al., 2015). A considerable amount of work has thus far been put into building and training statistical and machine learning models capable of accurate classifications of overall tissue samples as well as subsections of tissue samples (Amgad et al., 2019; Lu et al., 2020; Bian et al., 2021; Negahbani et al., 2021). In the context of cancer this work, often referred to as Histopathological Image Analysis or Digital Pathology, encompasses a broad collection of methods and goals ranging from developing models to assist in the quantitative scoring and staging of cancer biopsies to the classification of cells in tumor biopsies (Komura and Ishikawa, 2018). The value of these methods have important implications for basic science and translational research. Given the precision with which such algorithms can assess biopsies at a cellular level, it is now possible to uncover structures and associations which might not immediately apparent to human pathologists. This degree of precision also allows for rigorous modeling and examination of intratumoral heterogeneity with respect to their genomic determinants and clinical characteristics (Saltz et al., 2018; Li et al., 2019a).

Several research works in the past decade provide evidence regarding the importance of this line of investigation, especially for quantifying tumor-immune cellular interactions in cancer. Pagès et al. (2010) proposed immune reaction as the seventh hallmark of cancer, and laid out different associations between various types of immune cells and outcomes of interest. Among the different types of immune cells, the most well-studied are lymphocytes. Lymphocyte infiltration is both a meaningful prognostic indicator that can help inform treatment and predict survival across different types of cancer, including in colorectal cancer (Idos et al., 2020), breast cancer (Denkert et al., 2015; Dieci et al., 2018), and melanoma (Fu et al., 2019). This has led to increased interest in computational Tumor Infiltrating Lymphocyte (TIL) assessment, a sub-field of digital pathology devoted to developing computational methods to assess lymphocyte infiltration in biopsy images. Thus far, results that have emerged from this research have provided additional evidence regarding the importance of tumor infiltrating lymphocytes and their spatial features in the assessment of pathology images. Saltz et al. (2018) specifically examined the presence of lymphocytes across biopsies in several types of cancer, and found that certain summarization metrics of lymphocyte clusters were significantly associated with survival in certain types of cancer. Lu et al. (2020) also found that various spatial features of tumor infiltrating lymphocytes present in a sample of breast cancer biopsies were significantly associated with survival both marginally and after adjusting for other factors. Further, they found that these associations actually differed by tumor subtype. In addition to this, they found that certain spatial features of tumor infiltrating lymphocytes were also significantly associated with gene expression data.

From an analytical perspective, digital pathology research can be broadly classified along two axes: the type of quantitative models, and the scale of the assessment. Quantitative models range from primarily machine-learning in nature to primarily statistical. The most obvious way this difference manifests is in the models and methods used. Machine learning approaches tend to use tools like deep learning and convolutional neural networks and focus on predictive accuracy (Khosravi et al., 2018), while statistical approaches tend to use specialized model-based spatial methods and focus more on the quantification of uncertainty (Seal et al., 2022). The scale of assessment, on the other hand, relates to how much of the biopsy is used. Some methods operate on sub-regions of the biopsy images while others operate at the level of the entire biopsy (Ertosun and Rubin, 2015; Li et al., 2019b).

In theory, these two qualities-type of research and scale of assessment-are independent. In practice, methods that use machine learning approaches tend to operate on whole biopsy images or data, while statistical approaches tend to operate on sub-regions of images, or substantially smaller biopsies as with multiplex imaging data analysis (Seal and Ghosh, 2022). There has been limited work to build off of the cell classification abilities of modern machine learning methods at the level of full biopsy imaging. There are two probable reasons for this. First, many statistical methods designed for usage with spatial data rely on the definition of a window in which the analysis occurs. For biopsies, such a window is not readily defined, except in deliberately chosen sub-regions. While such issues can be avoided by using methods that do not rely on an explicit window, this tends to introduce additional computational challenges that may not scale to the biopsy level, especially as biopsy size and resolution increases. Second, even if one were to define a singular window for the entire biopsy, applying spatial methods across the entire space may prove problematic. In addition to being computationally intensive, it is not clear that measures of various local spatial characteristics summarized at the biopsy level are able to capture the heterogeneity across the biopsy, and thus may not be informative enough to be useful in practice. For these reasons, there has been limited work in development of spatial statistical models that use cell classification data from full biopsy images. To the author’s knowledge, the most notable attempts to do so have been Li et al. (2019b)’s work which utilized data classified by the ConvPath pipeline (Wang et al., 2019).

In order to address these issues, we propose SPARTIN (SPatial Analysis of paRtitioned Tumor-Immune imagiNg) pipeline. SPARTIN is unique in the spatial pathological imaging analysis literature in that it uses spatial statistical models to assess the association between tumor cells and immune cells across an entire partitioned biopsy, rather than select sub-regions. This allows for rigorous quantification of uncertainty in the style of Li et al., 2019b while still allowing for the assessment of full images in the style of Saltz et al. (2018). We accomplish this by partitioning the cell-level imaging data of each biopsy into non-overlapping sub-regions, which can then be modeled to capture the local infiltration patterns using techniques from spatial point process theory. Spatial point processes have long been used in the domain of ecology to rigorously investigate spatial relationships between various organisms (Högmander and Aila, 1999; King et al., 2012). More recently, methods from this field have been successfully applied within the biomedical domain to functional neuroimaging and magnetic resonance imaging data (Kang et al., 2011; Ray et al., 2015). Briefly SPARTIN uses whole slide images of biopsies that are partitioned into sub-regions using an iterative clustering algorithm and a Bayesian spatial marked point process model is then fit on each sub-region (see Figure 1). We discuss the constructions and justifications behind the choice of models underlying SPARTIN in Section 2. Using our model parameters, we construct a measure of local immune cell infiltration, termed Cell Type Interaction Probability (CTIP) that allows rigorous incorporation of uncertainty and is highly interpretable (on a probability scale) both within and across biopsies, and can be used to assess associations with genomic and clinical features. We evaluate the performance of SPARTIN on simulated data across multiple scenarios of cell interactions (Section 3.1) and show that our models are able to reliably distinguish positive interaction from negative and null interaction under a variety of scenarios (minimum AUC 0.85, maximum 0.97).

FIGURE 1

FIGURE 1. Overview of SPARTIN pipeline for single biopy image. Raw biopsy images (A) are preprocessed into cell types (tumor cells or immune cells, colored black and red respectively) and locations. The points are intensity thresholded, removing excess white space from the image and producing a more closely fitted window (B). Next, the biopsy is tessellated into “tiles” consisting of non-overlapping sub-regions containing comparable numbers of tumor cells (C). Bayesian Spatial point process models are fit on each tile (D), and tile-specific Cell Type Interaction Probability (CTIP) is calculated from the resulting posterior distributions (E), allowing us to capture heterogeneity in infiltration across the biopsy in a manner that can be aggregated (F) and appropriately summarized (G).

We demonstrate the utility of SPARTIN on a Skin Cutaneous Melanoma dataset consisting of 335 biopsies obtained from The Cancer Genome Atlas (TCGA) Genomic Data Commons Data Portal (Section 3.2). Our scientific goals are to quantify and characterize local spatial immune cell infiltration within and across melanoma biopsies and evaluate their association with genomic, phenotypic, and clinical outcomes. Using local spatial infiltration patterns quantified using CTIP, we found that CTIP was significantly negatively associated deconvolved immune cell prevalence scores including CD8+ T-Cells and Natural Killer cells. Furthermore, average CTIP scores differed significantly across previously established transcriptomic classes in melanoma (Akbani et al., 2015). CTIP was also more significantly associated with increased hazard of death, as compared to existing measures of spatial interaction. Beyond the inverse association with survival we found that CTIP was significantly inversely associated with clinical assessment of TIL abundance. Finally, we discuss the implications of this pipeline as well as potential future directions for the work (Section 4).

2 Materials and methods

2.1 Overview of SPARTIN

The overall schematic of the SPARTIN pipeline (for a single biopsy image) is shown in Figure 1. Briefly, for each whole-slide image of a biopsy (Panel A), the smallest window that fit all cells (Tumor and Immune) in the entire biopsy is constructed by thresholding the image intensity (Panel B). This window is then divided into “tiles” via a tessellation process that fully partitions the image into contiguous sub-regions containing similar numbers of tumor cells (Panel C; zoomed in version in Panel D). A Bayesian spatial point process model is then fit on each of these tiles separately, yielding a posterior distribution of the local interaction parameter (Panel E). For each tile, this local distribution is compared to a tile-specific null distribution to compute the local CTIP value, which is then combined to summarize the overall level of interaction at the biopsy level (Panel G). Each of these steps in detailed in the ensuing sections.

2.2 Intensity thresholding and partition

2.2.1 Motivation for partition

Each biopsy was partitioned into non-overlapping sub-regions, each of which was modeled separately as a point process. The primary motivation for structuring the analysis in this manner came from prior biological knowledge about tumor composition. Tumors are heterogeneous entities in many respects, and immune infiltration is no exception. Tumor biopsies can also be quite large relative to the scale on which interaction is defined. While we used a radius of interaction of 30 microns, some biopsies in our sample comprised tens if not hundreds of thousands of square microns of tissue. Additionally, setting aside practical computational issues with fitting a single model on an entire biopsy, such an approach would also both fail to capture spatial variability in immune infiltration and fail to allow for a more granular spatial investigation of the immune infiltration in a given biopsy. A natural way to improve over this method is to algorithmically partition the biopsy into non-overlapping sub-regions and fit models on each sub-region. Aside from better capturing the spatial variability of immune infiltration, an additional benefit of this method is that it allows for the parallelization of model fitting across the resulting sub-regions within a single biopsy. Much like selecting a quadrature when approximating an integral, calibrating the fineness of the partition entails a tradeoff between the precision with which one can assess the spatial heterogeneity within the biopsy and the computational load that a finer partition entails.

2.2.2 Partition pipeline

In order to partition a given biopsy into non-overlapping sub-regions, we began with the smallest bounding rectangular window that contained all cells. We then applied an intensity thresholding algorithm in order to find the smallest possible window that still contained all cells. This was accomplished by breaking the full image up into small (∼15 microns²) sub-windows, and computing a smoothed cell intensity on each sub-window. Sub-windows which had a sufficiently high intensity were combined to form the overall window on which the subsequent partition took place. The final window can thus be thought of as the union of many small square windows which, due to their size, are able to fit the contours of the differently shaped biopsies and exclude areas both outside and within the biopsies where no cells are present. The algorithm takes in three tuning parameters that determine the resulting partition. The same values for all three parameters were used across all biopsies, and the values themselves were selected through empirical tuning on a subset of the total dataset. For more details on the parameters themselves, see Supplementary Section S3.

Next, we applied a voronoi tessellation to the tumor cells within the intensity thresholded window, partitioning it into tumor cell specific sub-windows. We then applied a modified version of K-means to the tumor cells, such that each of the resulting clusters was constrained to be between a pre-defined range of cell counts. Finally, each tumor cell specific sub-window within a given cluster was combined into a tile, corresponding to each of the clusters from the K-means clustering; see Figure 1 for an illustration of the process. This results in a set of tiles that fully partition the intensity thresholded window. This partition uniquely defines the membership of each of the total cells (tumor or immune) into one of the resulting tiles. Because each resulting tile has a well-defined boundary and each cell in the biopsy belongs to exactly one tile, a hierarchical Strauss model can be fit on each tile to compute a tile specific value of each parameter in the model. Most notably, this allows for the estimation of a tile-specific CTIP value, which captures the local degree of tumor cell-immune cell interaction. By modeling this parameter locally to each tile, we are able to capture not only the overall level of infiltration in the biopsy, but also the potentially spatially heterogeneous nature of the infiltration.

It is worth emphasizing at this point the modular nature of our pipeline. The clustering method does not depend on the details of the model used, and the model in turn does not depend on details of the clustering method. Either can be exchanged for a different algorithm or model without disrupting the rest of the pipeline, so long as the output of the clustering algorithm is a spatial partition of the biopsy. Further, while the value ultimately selected to summarize the local infiltration at the tile level will obviously be informed by the exact model selected, there is flexibility in the choice of this quantity as well.

2.3 Spatial point process models

2.3.1 Data structure

Let N be the number of biopsies, each with N_i tiles. More formally, after the partition biopsy i can be represented by $T_{i1}\cup \cdots \cup T_{i{N}_i}$, where each T_ij closed subset of ${\mathbb{R}}^2$. Further, let x_ij = x_ij1 ∪ x_ij2, where each x_ijk is a set of points in $T_{ij}\subset {\mathbb{R}}^2$, and k = 1 for tumor cells and 2 for immune cells. Then for each element of the partition of each biopsy, this data can be naturally viewed as the realization of a marked point process. Broadly speaking, marked point process models are a family of models that treat the number of points, their locations, and their discrete classification (often called their “mark”) as random. We specifically chose to model this data as a hierarchical Multitype Strauss Process. Hierarchical Strauss Processes are unique in that they allow for the modeling of both the tendency of points of different types to be generally close (positive interaction) as well as generally far apart (negative interaction).

2.3.2 Model construction

The usage of a Strauss model to study spatial interaction in a marked point process is natural in this setting, since Strauss models were originally conceived for this purpose (Strauss, 1975). There are two primary model variants to consider when modeling multi-type data using the Strauss model: the hierarchical variant and the symmetric variant. These models are superficially similar, but differ in important ways. Broadly speaking the symmetric variant models the interaction between the different types of points jointly, whereas the hierarchical variant models the locations of certain types of points conditional upon the locations of other types. Using the hierarchical Strauss model thus requires an important choice for each pair of types as to which type of point will be “dominating,” i.e., conditioned upon, and thus unaffected by the other type. This choice, which must be made with a priori information about the nature of the process, has important consequences for the resulting model fit. Aside from the differing interpretations depending on the choice, different choices of dominating vs. non-dominating types will also result in different conditional likelihoods and thus parameter values (Högmander and Aila, 1999). It is also worth noting here that the term “hierarchical” is used differently in this context than it is often used in more typical statistical modeling. Rather than referring to mechanisms to account for dependence in observations or induce shrinkage, here the term refers to the hierarchical relationship between the dominating and non-dominating types: the locations of the dominating type do not depend on the locations of the non-dominating type, but the locations of the non-dominating type do explicitly depend on those of the dominating type.

We ultimately chose to use the hierarchical variant of the multi-type Strauss model over the symmetric variant for two key reasons. First, in our imaging application treating the locations of the immune cells as conditional upon the locations of the tumor cells is a priori biologically plausible. It has been established that there is an interplay between a tumor and the host’s immune system whereby the immune system responds to the tumor, thus shaping its development, and the tumor in turn “escapes” the immune system and influences its response (Khong and Restifo, 2002; Zhang and Zhang, 2020). While such dynamics are no doubt of interest, this natural evolution cannot be easily inferred from a “snapshot” in time that a biopsy presents. Despite this, conditioning on the locations of the tumor cells to model the immune response is a reasonable starting point: in order for the immune system to respond, there must at the very least be a tumor to respond to. The hierarchical Strauss Model thus serves as a first order approximation to this complex interplay. Second, the hierarchical Strauss model allows for the proper modeling of positive interaction between points of different types, while the symmetric model does not. Figure 2 depicts examples of different types of spatial interaction, ranging from negative to positive. The left and center panels depict the range of interactions that can be modeled by the traditional Strauss process: either the case where cells of both types actively avoid the other (left panel) or are spatially unaffected by the other’s presence (middle panel). Because we are interested in not just these cases but the case where points of both types tend to cluster around points of the other type (right panel), the hierarchical variant is the clear choice.

FIGURE 2

FIGURE 2. An illustration of spatial interaction between cells of different types. The leftmost panel depicts low/negative interaction; the middle panel depicts null interaction; the rightmost panel depicts high/positive interaction.

The density function of the hierarchical Multitype Strauss Model with two qualitative marks (tumor and immune cells in our setting) is defined by

$f\left({\mathit{x}}_1,{\mathit{x}}_2\right)\propto \exp \left(n_1{\beta }_1+n_2{\beta }_2+S_{R_{11}}\left({\mathit{x}}_1\right){\gamma }_{11}+S_{R_{22}}\left({\mathit{x}}_2\right){\gamma }_{22}+S_{R_{12}}\left({\mathit{x}}_1,{\mathit{x}}_2\right){\gamma }_{12}\right)(1)$

where in addition to x₁ and x₂ as defined above, n_t is the number of points of type t, β_t is the first order intensity of points of type t, $S_{R_{tl}}(\cdot )$ counts the number of pairs of points of types t and l within R_tl of one another where R_tl is selected a priori based on subject specific knowledge, and γ_tl captures the tendency of points of type l to be near points of type t. It bears mentioning that this model can be extended to arbitrary numbers of marks (i.e., types of cells). For our current application, we have focused on tumor-immune interactions, but it would be trivial to extend the chosen model to as multiple cell types as were available for a given data set.

2.3.3 Interpretation of parameters

As previously mentioned, γ_tl can be thought of as the degree to which points of type l tend to be close to or far away from points of type t. This allows us to distinguish between the mere relative abundance of different types of points (which is captured by the β₁ and β₂ parameters) and the actual spatial associations between the different types of points. This distinction is important in situations where there are substantial numbers of points of both types but no positive spatial association (and possibly a negative one). For examples, see Section 3.1 for results from the simulation study. When t = 1 and l = 2, this association is interpreted conditionally upon the locations of the type 1 points. Under the hierarchical Strauss model, γ₁₂ ∈ (−∞, ∞), with γ₁₂ ∈ (−∞, 0) implying negative interaction, γ₁₂ ∈ (0, ∞) implying positive interaction, and γ₁₂ = 0 implying no interaction at all. As one might intuitively expect, larger values of γ₁₂ correspond to more positive interaction, and smaller values correspond to a more negative interaction. Due to mathematical constraints on the density, we must have ${\gamma }_{11},{\gamma }_{22}\in \left(-\infty ,0\right]$. This constrains interaction between points of the same type to be modeled as negative or null, i.e., as in panels 1 and 2 of Figure 2. Since this constraint is unlikely to be satisfied in the context of a tumor biopsy, we set γ₁₁ = γ₂₂ = 0, in effect assuming no interaction between cells of the same type. Based on prior biological knowledge, we set R₁₂ = 30 microns. The need to select values of R₁₁ and R₂₂ is obviated by setting γ₁₁ = γ₂₂ = 0.

2.4 Model fitting

Except in the trivial case where γ₁₂ = 0, the normalizing constant of this distribution is computationally intractable. This motivates the usage of a pseudolikelihood formulation as outlined by Baddeley and Turner, 2000. For the density f(⋅) of the simplified hierarchical Strauss Model outlined in (1), define the conditional intensity function at a point u given θ = {β₁, β₂, γ₁₂} and points x observed in window A by

$\lambda \left(u|\theta ,\mathit{x}\right)=\left\{\begin{array}{cc}\frac{f\left(\mathit{x}\cup \left\{u\right\}\right)}{f\left(\mathit{x}\right)}\hfill & u\notin \mathit{x}\hfill \\ \frac{f\left(\mathit{x}\right)}{f\left(\mathit{x}-\left\{u\right\}\right)}\hfill & u\in \mathit{x}\hfill \end{array}\right.(2)$

Given this, the pseudolikelihood is defined as follows,

$PL\left(\theta |\mathit{x}\right)=\prod _{x_i\in \mathit{x}}\lambda \left(x_i|\theta ,\mathit{x}\right)\exp \left(-{\int }_A\lambda \left(u|\theta ,\mathit{x}\right)du\right)(3)$

In essence, this provides a computationally tractable alternative to the likelihood function that can be used in inference. Specifically, the integral in the pseudolikelihood function can be easily approximated by summing over a weighted quadrature on A, as described in Baddeley and Turner, 2000. See Supplementary Equation S1 for additional details.

Using this approximation of the pseudolikelihood function in place of the more standard likelihood function, Bayesian analysis can proceed in the style of King et al., 2012 by simply assigning priors to the parameters of interest and using techniques to sample from non-closed form posterior likelihoods. We assigned non-informative normal priors with mean 0 and variance 10⁶ to γ₁₂, β₁, and β₂. In all analysis presented below the quadrature and weights used to estimate the integral in the pseudolikelihood function was generated by the spatstat package (Baddeley and Turner, 2005). Samples from the posterior were taken using JAGS via the R2jags package (Yu-Sung et al., 2015); see Supplementary Algorithm S1 for details.

2.5 Cell type interaction probability

2.5.1 Motivation

For ease of exposition, we suppress the subscript in γ₁₂ parameter, and refer to it as γ_• in the following sections. γ_•, β₁, and β₂ allow us to distinguish between the tendency of cells of different types to be spatially near one another (captured by γ_•) and the relative abundance of the cells of each type (captured by β₁ and β₂), as well as the uncertainty around these various tendencies. However, this does not capture the difference between the observed spatial association and what one would expect to observe by chance given a particular configuration of tumor cells. Even if the locations of the immune cells are completely random and unrelated to those of the tumor cells, there is still a possibility of a configuration that may indicate positive interaction when considered in isolation. In order to properly account for this, the observed distribution of γ_• must be compared to a counterfactual distribution that captures the tendency when there is no interaction. This motivated the development of Cell Type Interaction Probability (CTIP).

2.5.2 Definition

Let π(γ): = f(γ_•|x₁, x₂) be the posterior distribution of γ_• conditional upon the observed tumor cells and immune cells, and ${\pi }_0(\gamma ):=f({\gamma }_0|{\mathit{x}}_1,{\tilde{\mathit{x}}}_2)$ be the posterior distribution of the interaction parameter conditional upon the observed locations of the tumor cells and ${\tilde{\mathit{x}}}_2$, a set of points distributed as a Poisson Process. Further, assume independence between the true distribution and the null distribution. Then we define the CTIP (denoted by r) as:

$r=\iint I\left(\gamma >{\gamma }_0\right)\pi \left(\gamma \right){\pi }_0\left({\gamma }_0\right)d\gamma d{\gamma }_0,(4)$

where I(⋅) denotes the indicator function, and both π(γ) and π₀(γ₀) are as defined above. By construction, the observed parameter is independent from the null parameter, since immune cells drawn from a Poisson process should not yield any information about immune cells realized under either positive or negative interaction. More concretely, by the definitional independence of ${\tilde{\mathit{x}}}_2$ from x₂, it follows trivially that $f({\gamma }_{•},{\gamma }_0\mid {\mathit{x}}_1,{\tilde{\mathit{x}}}_2,{\mathit{x}}_2)=f({\gamma }_{•}\mid {\mathit{x}}_1,{\mathit{x}}_2)f({\gamma }_0\mid {\mathit{x}}_1,{\tilde{\mathit{x}}}_2)$. This construction means that CTIP is the expectation of an indicator function, and thus can be interpreted as a probability. Specifically, it represents the probability that the observed value of γ_• is larger than what we would observe if there was truly no underlying interaction.

2.5.3 Sampling from the empirical null distribution

The problem of estimating CTIP hinges on being able to sample from the distribution of γ₀, since estimation of the posterior can proceed as described in Section 2.3. In order to estimate the null distribution of γ₀ conditional upon the location of the observed tumor cells, we used simulation to construct an empirical null distribution. Recall that when there is no interaction, by definition γ = 0. This means that the hierarchical Strauss given in Eq. 1 reduces to two independent Poisson processes with intensities β₁ and β₂, which are trivial to simulate.

Given this, the estimation of the null distribution for a given tile proceeded as follows. First, the first order intensity of the immune cells was estimated using the standard estimator, ${\widehat{\beta }}_2=\log (\frac{I_{ik}}{A_{ik}})$ where I_ik is the number of immune cells observed in tile k of biopsy i and A_ik is the area in square microns of the tile. Second, S simulations of immune cells were generated from a poisson process with intensity ${\widehat{\beta }}_2$, S being selected a priori. Third and finally, each simulated set of immune cells was superimposed over the actual observed tumor cells, and samples were drawn from the resulting posterior distribution. These samples, combined across simulations 1, …, S approximate π₀(γ) as previously defined.

2.5.4 Computation of CTIP

Having established the ability to sample from the posterior distribution as well as the null distribution for γ, estimation of CTIP can be done via stochastic integration. Note that Eq. 4 is equivalent to $E_{\gamma ,{\gamma }_0}[I(\gamma >{\gamma }_0)]$. Thus, given P posterior samples γ₁, …, γ_P from the posterior π and γ₀₁, …, γ_0P from the empirical null π₀, (5) can be estimated by

$\widehat{r}=\frac{1}{P}\sum _{j=1}^{P}I\left({\gamma }_j>{\gamma }_{0j}\right)(5)$

2.5.5 Interpretation

By definition, r ∈ [0, 1]. As previously mentioned, r can be construed as the expected value of an indicator random variable with respect to the joint distribution of π(γ) and π₀(γ₀). This naturally yields an interpretation of r as a probability. Specifically, it can be thought of as measuring the posterior probability that the observed value of γ_• is larger than the value that would be observed given the observed tumor cells and immune cell intensity if there were no underlying interaction.

2.5.6 Model outputs and summaries

Because models can be fit for each sub-region T_ij of biopsy B_i, CTIP can also be computed for each tile. The result for each biopsy is a collection of tile specific estimations of CTIP for each biopsy, $\{{\widehat{r}}_{i1},\dots ,{\widehat{r}}_{i{N}_i}\}$ for i = 1, …, N. These estimates are then summarized at the biopsy level into a singular value that quantifies the overall level of infiltration in the biopsy, which can be used in downstream analyses to investigate the association between biopsy level of tumor-immune interaction and other patient features of interest-see Section 3.2 for results. We chose to use the mean value, since even under the assumption of spatial association between tiles the result is an unbiased estimate of the biopsy level mean. Because the values themselves encode the local spatial association for each tile, the result is an appropriate proxy for the amount of interaction between tumor cells and immune cells at the biopsy level. We further logit-transformed the tile level CTIP values prior to computing the biopsy level mean in order to make the resulting distribution across biopsies closer to a Gaussian distribution. See Algorithm 1 for the full SPARTIN method.

Algorithm 1. SPARTIN algorithm.

Input: N biopsies

for biopsy i = 1, …, N do

Intensity threshold biopsy i

Partition intensity thresholded biopsy i into N_i tiles

for tile j = 1, …, N_i do

Estimate CTIP ${\widehat{r}}_{ij}$ for tile j

end for

Compute ${\widehat{m}}_i=\text{mean}(\text{logit}({\widehat{r}}_{i1}),\dots ,\text{logit}({\widehat{r}}_{i{N}_i}))$

end for

return ${\widehat{m}}_1,\dots ,{\widehat{m}}_N$ for usage in downstream analysis

3 Results

3.1 Simulations

3.1.1 Simulation overview

To compare our model’s detection of spatial association between different cell types, we conducted simulation studies in which we tested the ability of SPARTIN to accurately identify positive interaction across a range of different simulated cell compositions and degrees of spatial association. For each of four sets of simulations, the number of simulated tumor cells and immune cells were set a priori at T_s and L_s. These levels across the four simulations were, respectively, 20 and 50, 50 and 20, 50 and 50, and 100 and 100. Further, interaction was controlled by a parameter ϕ ∈ [−1, 1], with −1 indicating the most negative possible interaction using our simulation method and 1 indicating the most positive possible interaction using our simulation method. For each of the four settings, 50 data sets were simulated for ϕ = ±1, ±0.8, ±0.6, ±0.4, ±0.2, and 100 were simulated for ϕ = 0. Because the goal was classification as either positive interaction or non-positive (i.e., null or negative) interaction, we used AUC as our summary metric. This is a natural choice, since CTIP is interpreted as the posterior probability of positive interaction. As a benchmark for performance, we compared CTIP to the G-cross function, a commonly used non-parametric statistical method for assessing the spatial relationship between different types of points in marked point processes (Van Lieshout and Baddeley, 1999).

3.1.2 Simulations of negative interaction

For a given combination of T_s and L_s, the simulations of positive interaction and simulations of negative interaction proceeded differently. For the simulations of negative interaction (ϕ ∈ [−1, 0]), T_s tumor cells and L_s immune cells were simulated as independent poisson processes in regions of the window that overlapped to varying degrees. The overlap was controlled by ϕ, such that the processes overlapped on (100 ⋅ (1 + ϕ))% of the window in which the simulation occurred. Note that when ϕ = −1 there was no overlap, and when ϕ = 0 there was complete overlap and the simulation reduced to generating two independent Poisson Processes within the same window.

3.1.3 Simulations of positive interaction

For the simulations of positive interaction ($\varphi \in \left(\mathrm{0,1}\right]$), T_s tumor cells were simulated under a Poisson process. After their locations were determined, L_s immune cells were simulated. For each immune cell l_i, a Bernoulli random variable C_i∼Bern(ϕ) was drawn. If C_i = 1, l_i was simulated within 30 microns of a randomly selected tumor cell. Otherwise, l_i was simulated from a Poisson process. Thus, the level of interaction was again controlled by ϕ, with ϕ = 0 again corresponding to two independent Poisson Processes and ϕ = 1 corresponding to a situation in which all immune cells are within 30 microns of at least one tumor cell.

3.1.4 Simulation assessment

Because CTIP is most naturally interpreted as a probability, in order to classify a given collection of points as indicating positive or negative interaction in practice a cutoff value must be chosen above which CTIP is considered to indicate positive interaction and below which it is not. This leads to a natural tradeoff: if the cutoff is too low, then many cases where there is no positive interaction will be incorrectly classified as having positive interaction. If the cutoff is too high, the opposite will occur. This tradeoff can be visualized across different thresholds using a receiver operating characteristic (ROC) curve, which can in turn be summarized by the area under the curve (AUC). AUC ranges between 0.5 and 1, with 0.5 indicating total inability to correctly discriminate between different outcomes and 1 indicating perfect ability to discriminate. AUC is therefore a natural choice for a metric to assess the ability of CTIP to distinguish between positive and non-positive interaction between cells of different types.

As a benchmark for performance, we examined the performance of the G-cross function in discriminating positive interaction from non-positive interaction. The G-cross function is a non-parametric estimate of the cumulative distribution functions of minimum distances from cells of one type to those of another. In this application, the G-cross function evaluated at a specific distance yields a non-parametric estimate of the probability that a tumor cell will have at least one immune cell within that distance. We chose to evaluate the G-cross function at a distance of 30 microns in order to be consistent with the modeling of CTIP.

3.1.5 Results

For examples of simulated data as well as visualizations of results for CTIP, see Figure 3. Across the different simulation settings, accurate classification was achieved using CTIP. The top panels depict a sample of simulations across a range of interaction values, from ϕ = −1 to ϕ = 1. The bottom left panel shows the distribution of the estimated CTIP values for each set of simulations across each of the four sets of simulations and the values of ϕ, plotted on the logit scale. The bottom right panel shows the AUC for each of the four simulation settings across all values of ϕ. The results of these simulations demonstrate that CTIP reliably distinguishes between a broad ranges of interaction levels, across a similarly broad range of relative abundances in cell types. The lowest overall AUC (0.85) was achieved in simulation 3, while the highest was achieved in simulation 1 (0.96). For the G-cross function, the highest overall AUC was achieved in simulation 2 (0.8), while the lowest was achieved in simulation 1 (0.73)—see Supplementary Figure S4 in the supplementary materials for more detailed results for the G-cross function, and Supplementary Section S5 for additional simulation results examining the relative performance when model assumptions are violated.

FIGURE 3

FIGURE 3. Simulated examples and results. Top row: examples of simulations of positive and negative interaction for different values of ϕ, all from simulation setting 1. Bottom left: box plots of CTIP across all simulated compositions for different values of ϕ on logit scale. Bottom right: ROC curves and associated AUC for classification (positive interaction vs. non-positive interaction) for different simulation settings and the two methods investigated (CTIP and G-cross).

3.2 Application

We applied the SPARTIN pipeline to a data set consisting of 335 high definition images of SKin Cutaneous Melanoma (SKCM) biopsies stained using hematoxylin and eosin, obtained from The Cancer Genome Atlas (TCGA) Genomic Data Commons Data Portal. SKCM is an appealing target for the investigation of immune cell infiltration for several reasons. SKCM has been shown to be particularly responsive to immunotherapy in some cases (Achkar and Tarhini, 2017; Franklin et al., 2017). It is possible that the ability to quantify immune infiltration at a large scale may allow for more detailed investigation into the scenarios in which this treatment may be most effectively deployed. It is also well established that SKCM has an unusually high mutational load amongst the various cancer types (Berger et al., 2012). The ability to quantify spatial infiltration patterns may allow for further downstream investigations into not only genomic associations of this occurrence but associations with other phenotypic and clinical outcomes as well.

3.2.1 Cell classification model

3.2.1.1 Whole-slide image (WSI) retrieval and obtaining training labels

WSI from 20 different patients were used as our primary training and testing data. The images were imported from The Cancer Genome Atlas (TCGA) Genomic Data Commons Data Portal, which is a repository of validated datasets from various National Cancer Institute Programs: https://www.cancer.gov/tcga. Also known as the gigapixel pathology image, one image slide contained two tissue smears stained with hematoxylin and eosin, imported in a vendor specific format. The imported image for classification is scanned at the “high” magnification level in a microscope of ×40. The pixel size for images was approximately 0.25 microns per pixel. Manual marking of a minimum of ten samples from each of the classes, i.e., Tumors, Immune, Macrophage, and other cell was performed by a pathologist on all slides. Care was taken to ensure that the samples were as diverse as possible, to account for all possible morphologies of the same cell type across all cases. In total, 1,250 annotations were made in all the images, with each annotated cell labelled by a pathologist. The reasons for the comparatively smaller labelled dataset available lie in the large size of the slide, and pathologist availability; the large size makes parsing through difficult and explains the diversity of labelled structures in each image.

3.2.1.2 Whole-slide image (WSI) nucleus segmentation

Before nuclear segmentation and extraction can be performed on the WSIs, certain pre-processing steps are carried out to ensure staining uniformity across images. A representative and generalizable stain vector estimation for each of the stains being used in the image is estimated, to normalize the staining intensities across all the images. In our application, vectors from hematoxylin, eosin and residual stains were estimated from a standard image, selected by the clinician. These vectors are then applied across all images to keep the stain detection parameters uniform across the dataset (Bankhead et al., 2018).

Due to the variations in staining procedures across the dataset, it is difficult to accurately isolate whole cells on the slide, as the cytoplasmic boundaries are not well-defined due to lack of membrane staining. As heterogeneity in nuclear morphology has been shown to be a good discriminator, only the nucleus is segmented out of each cell in every image for classification purposes, using Qupath (Yuan et al., 2012). Watershed segmentation is used to obtain separated and contoured nuclei from the hematoxylin color channel of the whole slide image, in a patch-based manner (Bankhead et al., 2018). The separate hematoxylin channel is obtained by performing color deconvolution to separate out the stains used in the slide (Ruifrok and Johnston, 2001). It was observed that an average of 200,000 cells were segmented out per image, and morphological and intensity features such as spatial location, eccentricity, circularity, and stain intensities were computed. All pertinent image analysis and nuclear segmentation was performed using Qupath, an open-sourced software platform that can be used for a range of pathological image analysis applications (Bankhead et al., 2017).

A Random Forest model is a type of ensemble learning method, where the weak predictive power of multiple decision trees is aggregated to produce an accurate result (Breiman, 2001). Our four-class classifier model was developed with the pathologist labelled dataset of 1,250 morphologically diverse cells, each belonging to Tumor, Immune, Macrophage, or other cells (including cell types such as stromal cells). 5-fold cross-validation was used to assess model accuracy and the receiver operator characteristic curve AUC for each class, with proportional representation of all 4 classes ensured. After adjusting for multiple parameters, including the number of decision trees, the accuracy was obtained in the range of 87%–91%.

This step allowed for a preparation of a dataset for each cell with its spatial location using global coordinates and class of the cell, which is then used for further downstream analysis. The training, testing, and classification of cells from the WSIs was performed on MATLAB version 2017A. For visualizations of distributions of cell counts across biopsies, see Supplementary Figures S1–S3. For this analysis, only data for tumor cells and immune cells were kept to quantify the immune cell infiltration using SPARTIN.

3.2.2 SPARTIN analyses of SKCM

While we are broadly grouping immune cells together for the purposes of this analysis, it is important to note that within the microevironment of melanoma there is a great deal of heterogeneity in the types of immune cells present. Tumor Infiltrating Lymphocytes (TILs) receive a disproportionate amount of attention in cancer literature due to the well-established prognostic value of the presence. However, in the context of melanoma this focus can obscure the tremendous amount of variability in the types of immune cells present in the microenvironment. Even referring to TILs as a monolith obfuscates the different cell types captured under this umbrella, many of which serve different purposes (Mukherji, 2013; Weiss, 2016; Antohe et al., 2019). The richness of the immunological aspect of the melanoma microenvironment makes it a logical target of analysis via the SPARTIN pipeline.

3.2.2.1 CTIP quantification

In order to compute CTIP, models were fit in using the methods described in Section 2. The same settings were used for all biopsies across the pipeline. Each biopsy was partitioned such that the resulting tiles would have approximately 75 tumor cells per tile, with the minimum possible tumor cells per tile being 50 and the maximum being 100. For an analysis of the sensitivity to the biopsy level estimation of CTIP to the number of tumor cells per tile as well as the clustering algorithm used, see Supplementary Section S6. The CTIP quantifications are shown in Figure 4 Panel A wherein we have highlighted three biopsy images (top row). In the middle row, the color of each tile indicates the value of CTIP estimated for that tile. More varied colors across a given biopsy are indicative of more spatial variation in CTIP, and thus infiltration. CTIP was summarized at the biopsy level by taking the empirical mean across all tiles for a given biopsy, which are plotted in the bottom row. As can be seen, the range of CTIP was 0.26–0.99 with a median CTIP value of 0.69 and an interquartile range of 0.19. This points to the fact that SKCM has high level of immune cell infiltration and these results are consistent with the conventional wisdom that melanoma is a generally more immunogenic cancer relative to other variants (Mukherji, 2013; Fu et al., 2019).

FIGURE 4

FIGURE 4. Visualization of selected biopsies, and distribution of CTIP across all biopsies analyzed. Selected examples of actual biopsies, both raw spatial data and with tessellations and CTIP by tile. The leftmost biopsy pictured has among the lowest mean estimated CTIP in the sample (0.4), while the rightmost has among the highest (0.94). The middle biopsy has average mean CTIP (0.7), but also exhibits more spatial variation than the left or right biopsy. Bottom: distribution of mean CTIP by biopsy for all biopsies investigated. This illustrates the between-biopsy variation in overall levels of infiltration.

3.2.3 Genomic associations

We conducted a focused analysis on CTIP and genes that are associated with immune activity. Bhattacharya et al., 2018 have collected and classified a list of 1,793 unique genes associated with various aspects of human immune activity. Notably, this includes genes related to processes that involve or are directly related to immune cells of interest including CD8⁺ T-cells, Natural Killer cells, and B cells. Of these, gene expression data was available for 1,305 genes across 330 patients. We assessed the univariate association between biopsy level mean logit-CTIP and these genes using Spearman correlation. The advantage of Spearman correlation as opposed to the more standard Pearson correlation is that the former does not assume a linear relationship between the underlying variables of interest. Such assumptions can be problematic, particularly when there is no strong reason a priori to believe the relationship between the two variables is of a particular form. However, like Pearson correlation Spearman correlation is defined to lie in [−1, 1], with each extreme indicating the same directionality and strength of association as Pearson Correlation. After applying a Bonferroni correction, we found that 28 genes were significantly associated with CTIP. See Figure 5 for a volcano plot of results; for the complete list of genes, see Table 1. Most notably, mean logit-CTIP was most significantly negatively associated with CD244, a gene associated with increased cytotoxicity in natural killer cells (Agresta et al., 2018). In addition, mean logit-CTIP was significantly negatively associated with expression of Cathepsin S, a gene associated with immunosuppression in the tumor microenvironment, and C-C motif chemokine receptor 8 (CCR8), a gene related to immunosuppressive T regulatory cells (Fuchs et al., 2020; Weaver et al., 2022).

FIGURE 5

FIGURE 5. Volcano Plot of Spearman correlations between CTIP and gene expression values. Negative log₁₀ p-values plotted against Spearman Correlation for ImmPort genes. Colors indicate ImmPort Gene List, which are the collections provided by ImmPort to organize genes into groups of similar function; Red line indicates Bonferroni correction-see Table 1 for full list of significant genes.

TABLE 1

TABLE 1. Full ImmPort Analysis Results. ImmPort Genes significantly associated with average logit biopsy level CTIP after Bonferroni correction, alpha = 0.05.

We also investigated the association between CTIP values and tumor deconvolution data in order to assess the association between the presence of specific immune cell types and overall tumor-immune interaction. Because the specific immune cell types were not ascertained in the data collection, deconvolution algorithms can offer an estimation of the relative prevalence of different types of immune cells. We obtained data from the Tumor Immune Estimation Resource 2.0 (TIMER 2.0): http://timer.cistrome.org/. Using data from TIMER2.0, we examined the association between the prevalence of different types of immune cells including CD8⁺ T cells, Natural Killer cells, and B cells with biopsy level mean logit-CTIP. Among the available options from TIMER 2.0, we decided to use the MCP-counter algorithm (Becht et al., 2016) based on the analysis of Sturm et al., 2019, since it was judged to be most effective in detecting the presence and prevalence of the most relevant types of immune cells. We investigated the association of the score of each type of immune cell estimated by MCP-counter with biopsy level mean logit-CTIP using Spearman correlation (as done previously for gene expression). Significance was assessed using the standard test of statistical significance of Spearman correlation.

After applying a Bonferroni correction (α = 0.05), we found that six different immune cell scores as computed by MCP-counter were significantly negatively associated with biopsy level mean logit-CTIP: CD8⁺ T cells, B cells, Monocytes, Macrophages, Myeloid Dendritic Cells, and Natural Killer cells. No cell types were significantly positively associated with biopsy level mean logit-CTIP after the Bonferroni correction, though the magnitude of the positive association with Cancer Associated Fibroblasts (CAFs) is notable, and while not statistically significant still highly consistent with a truly positive underlying association between biopsy level CTIP and prevalence of CAFs. See Figure 6 and Table 2 for full details on specific cell type score associations. CD8⁺ T cells, B cells, and Natural Killer cells are generally associated with better prognosis in SKCM. This would imply, perhaps somewhat counterintuitively, that higher mean logit-CTIP would be expected to be associated with poorer survival, and lower levels of pathology assessment of lymphocyte infiltration. However, this is precisely what we observed, as discussed in next Section.

FIGURE 6

FIGURE 6. Volcano plot of association with deconvolution data. Negative log-10 p-values plotted against Spearman Correlation between CTIP and cell type prevalance scores estimated by MCP. Red line indicates Bonferroni correction, alpha = 0.05.

TABLE 2

TABLE 2. Associations with cell prevalence. Cell type associated with average logit biopsy level CTIP. p-values are presented uncorrected; see Figure 6.

3.2.4 Association with phenotypic and clinical outcomes

3.2.4.1 Transcriptomic classes and pathology asessments

Akbani et al. (2015) also identified three transcriptomic classes by applying consensus hierarchical clustering techniques to gene expression data from 1,500 genes: the “immune” subclass, the “keratin” subclass, and the “MITF-low” subclass. Most notably for our current application, the immune subclass was characterized by overexpression of genes associated with T cells, B cells, and Natural Killer cells. Of the patients classified using these methods, 235 were present in our data set. Of these 235, 114 (49%) were in the immune subclass, 78 (33%) were in the keratin subclass, and 43 (18%) were in the MITF-low subclass. Figure 7 demonstrates a clear visual difference between the immune subclass and both the MITF-low and keratin subclasses. In order to determine if these apparent differences were statistically significant, we analyzed the differences in average logit-CTIP between each pair of classes using a standard two-sided t-test, and found that the mean logit-CTIP in the immune class was significantly lower than either the keratin or the MITF-low subclass (p < 10^–5 and p = 0.011, respectively). We found no significant difference between the mean logit-CTIP in the MITF-low subclass and the keratin subclass (p = 0.21). This is highly consistent with the results of the deconvolution analysis, which suggest that the degree of abundance of such cells are inversely associated with mean logit-CTIP.

FIGURE 7

FIGURE 7. Average CTIP by transcriptomic subclass and overall survival hazard ratios. (A) Average logit-CTIP by transcriptomic class according to Akbani et al. (2015). (B) Hazard ratios from Cox Proportional Hazards models. Points at 1 indicate reference category. Age was standardized by subtracting the mean value and dividing by the standard deviation. Logit-CTIP was significantly associated with increased hazard of death (p = 0.019), while normalized logit-MCON was non-significantly associated with an increased hazard of death (p = 0.057) and normalized logit-G-cross was non-significantly associated with a decreased hazard of death (p = 0.32).

In addition to classifying the biopsies into transcriptomic subclasses, Akbani et al., 2015 had pathologists assess biopsies for lymphocyte infiltration. Biopsies were scored from 0 to 3 on lymphocyte distribution, with 0 indicating no lymphocytes present in the tissue and 3 indicating that lymphocytes were present in over 50% of the tissue. They were also scored from 0 to 3 on lymphocyte density, with 0 indicating an absence of lymphocytes and 3 indicating a “severe” presence. These measures were added to create a Lymphocyte Score, a measure that ranged from 0 to 6 meant to summarize the general presence and degree of lymphocyte infiltration in that biopsy. Of the 235 patients in our sample that were assessed by the pathologists, 29% had a score of zero, 23% had a score of two, 10% had a score of three, 13% had a score of four, 15% had a score of five, and 9% had a score of six. Note that by definition of the component scores, a score of one is not possible. We performed a linear regression of Lymphocyte Score on mean biopsy level logit-CTIP, treating Lymphocyte Score as continuous. We found that a one unit increase in mean biopsy level logit-CTIP was highly significantly associated with a 0.45 unit decrease in mean Lymphocyte Score (p ≪ 0.0001). This result is also consistent with both previous results, and provide further evidence for the apparent negative association between mean logit-CTIP and the presence of specific immune cell types generally associated with positive prognosis.

3.2.4.2 Survival analysis

In order to assess association between logit-CTIP and overall survival, we fit a Cox Proportional Hazards model using clinical data as well as CTIP values. Clinical patient data were retrieved from the TCGA website and matched to biopsy images via TCGA identifier. In addition to adjusting for mean logit-CTIP, we adjusted for cancer stage, age in years, and sex. Age in years was standardized by subtracting the mean value and dividing by the standard deviation.

Using this model, we found that after adjusting for other factors an increase in normalized logit-CTIP was significantly associated with increased hazard of death (p = 0.02). See Figure 7 for hazard ratio estimates and associated confidence intervals. In addition to this, we found that normalized age was significantly positively associated with an increased hazard of death (p < 0.05), as was having stage 3 disease (relative to stage 1, p < 0.05). Having stage 4 disease was positively associated with increased hazard of death, though not significantly. This is most likely due to the lack of patients with stage 4 disease in the data set.

In order to assess the performance of CTIP relative to a more traditional measure of spatial association the same model was fitted with average logit-CTIP exchanged for the average value of the estimated G-cross function (discussed previously in Section 3.1) as well as the Mark Connection Function (MCON). The Mark Connection Function is another natural point of comparison because it is a measure that is commonly used to investigate spatial associations between points with discrete marks, it is constrained to lie in [0,1] like CTIP, and it has been used elsewhere as a comparison point for similar survival modeling using spatial information (Li et al., 2019b). To produce a single value, the MCON function was evaluated at r = 30, which is the same as the radius of interaction used in our computation of CTIP as well as the G-cross function. Both G-cross and MCON values were logit-transformed. Because the values remained skewed even after logit-transform, they were subsequently centered and scaled by subtracting their respective mean and dividing by their respective standard deviation. While increased normalized logit-MCON was also associated with an increased hazard of death, the association was not significant (p = 0.057). The G-cross function was associated with a decreased hazard of death, though the association was also not significant (p = 0.32). Figure 7 summarizes the results across all three models.

This result, along with the other results presented in this paper, provide evidence of an immune phenotype in SKCM that corresponds to poor prognosis from a number of different perspectives. With respect to overall survival, we have demonstrated that adjusting for other relevant clinical factors, an increase in CTIP is associated with an increased hazard of death. For gene expression data, average CTIP is associated with decreased expression of genes related to immune cells that are generally associated with good overall prognosis. And finally, we found that average CTIP is significantly lower in biopsies that have been classified as immune enriched through gene expression clustering, as well as biopsies that have been assessed by pathologists to have higher densities and spatial distributions of tumor infiltrating lymphocytes. These results fundamentally support the notion that CTIP is capturing an anti-immunogenic phenotype that is associated with poor prognosis, possibly due to an overall negative association with the presence of immune cells that are associated with improved prognosis such as CD8⁺ T Cells and Natural Killer cells.

4 Discussion

The SPARTIN method provides a general framework for applying traditional spatial statistical techniques to whole biopsy imaging data. More specifically, it provides a modular algorithmic mechanism for partitioning the biopsies into more readily analyzable sub-regions that can be assessed for local spatial patterns using marked point process models. We have demonstrated the utility of such a method using a novel characterization of spatial cellular interaction, CTIP, that allows rigorous quantification of uncertainty and is highly intepretable. We demonstrate through simulation that CTIP can reliably distinguish between different patterns of spatial interaction across a range of different dependencies. We also demonstrate the utility of SPARTIN through a comprehensive analyses of an SKCM dataset, wherein we quantify and evaluate the spatial extent of tumor-immune infiltration and its association with genomic, phenotypic and clinical outcomes. These results are particularly suggestive in the context of SKCM. Melanoma has a uniquely rich literature regarding the potential applications of immunotherapy. While there are many factors that have contributed to this, not the least of which is the relative accessibility of skin cancers compared to other types of cancer, the net result is a far more mature understanding of immunotherapy in the context of melanoma than other types of cancer (Oble et al., 2009; Nguyen et al., 2010; Pitcovski et al., 2017). Thus, possible implications of the spatial immune organization in melanoma for immunotherapy may be more readily interpreted within the context of Melanoma than other cancer types, particularly as the ability to distinguish between different types of immune cells improves.

The SPARTIN pipeline represents a valuable contribution in and of itself to digital pathology through its ability to model and quantify immune infiltration across entire biopsies in a way that captures meaningful variation along with uncertainty quantifications across patients. However, it also creates numerous opportunities for future work in this area. As algorithms for cell-level image classification improve, the opportunities for more and more detailed quantitative analysis of histopathological imaging data will become both more numerous and more fruitful. Specifically, as the ability to reliably distinguish between different types of immune cells will allow for more granular investigations into the prognoses associated with interaction between immune cells of different kinds and the tumor. Moreover, as spatial biology techniques (e.g., spatial multiplex imaging and spatial transcriptomics) evolve and data becomes abundant, so too will opportunities for synthesizing data on the relative spatial locations of different cell types along with cell or spot specific gene/protein expression data through rigorous, principled, and complex modeling. It would be straightforward (particularly in a Bayesian framework) to model a function of one or all of the parameters of interest as a linear combination of the local genes/proteins that are known to be relevant to immune response. Alternatively, future methods in this area could be developed at the intersection of spatial data and high-dimensional omic features.

The SPARTIN pipeline could also be generalized for further investigation not directly related to assessing immune infiltration. Partitioning each biopsy into non-overlapping sub-regions invites the application of other tools from the spatial statistical catalogue to other relevant spatial features of the tumor microenvironment. In fact, SPARTIN can be utilized even when investigating features of the tumor that are not spatial in nature. So long as a feature (such as average tumor cell size) can be quantified at the cellular level and summarized locally at the level of sub-regions, SPARTIN provides a framework for mapping the variation of that feature across the entire biopsy.

One limitation of this analysis is that immune cells were modeled as a single agglomerated class rather than as separate cell types belonging to the same family. While our results paint a consistent picture from several different angles, the exact mechanism of these various associations could be further refined by utilizing a more precise categorization of the various types of immune cells. It is worth emphasizing, however, that nothing about the structure of our pipeline depends on the presence of only one type of immune cell in addition to the cancer cells. The SPARTIN method can be generalized to multiple types cells through multicategory marked point process models, though this would raise additional computational and practical challenges. Another potential limitation of our analysis that suggests an avenue for future methodological development is the presence of repeated measurements within biopsies. Note that the result of fitting models on the partitioned tumor microenvironment is a set of repeated observations with some underlying spatial covariance structure related to the spatial relationships between the different subspaces. For now, we have ignored this spatial covariance, and instead chosen to summarize the variation at the biopsy level. While this does take into account local spatial dependencies (i.e., within tiles), it is possible that incorporating additional global spatial structure across the biopsy (i.e., between tiles) could lead to additional gains in efficiency, and possibly insight into the nature of the spatial variation in immune activation in the tumor microenvironment. Thus, future work will likely include the development of methods to efficiently analyze this complex structured areal data. Finally, it should be noted that the hierarchical Strauss model makes two important assumptions about the nature of the data, those of stationarity and isotropy. This essentially amounts to the assumption that on the level of the tiles for each biopsy on which the models are fit, the point process would be unaffected by horizontal translation (stationarity) or rotation (isotropy). Because tumor biopsies do not have a natural orientation in ${\mathbb{R}}^2$, it is plausible to assume isotropy. Stationarity is somewhat more challenging to justify a priori. However, because the model fitting takes place on relatively small sub-sections of the biopsy, i.e., the “tile” specific level, it is quite plausible that this assumption is at least locally satisfied well enough. Since each tile is modeled using separate spatial parameters, this accounts for non-stationarity at a biopsy level and the ultimate measure of interest, CTIP, still yields a reliable estimate of the degree of interaction between cell types. Still, future work could involve employing methods that relax or eliminate these assumptions during the modeling of interaction.

Software is available at https://github.com/bayesrx/SPARTIN. The website for the project, which features an application that can be used as a visualization companion to the software can be found at https://nateosher.github.io/SPARTIN. Using the SPARTIN R package, data can be exported and visualized on the web application in order for users to interactively navigate the patterns of interaction at both the biopsy and individual tile level.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Author contributions

NO and VB conceptualized the project. NO developed the idea for the partitioning of the biopsy and for cell type interaction probability, and wrote the introduction, methods, results, and conclusion section. VB and JK oversaw the creation, implementation and visualization of the analyses pipeline. AR oversaw the preprocessing and preparation of the data by SK, and SK also contributed to writing the results section. AR and SK also provided guidance on the interpretation of results. VB and JK reviewed the manuscript and provided feedback for its ultimate construction and form.

Funding

VB was supported by NIH grants R01CA244845-01A1, and P30 CA46592 (Cancer Data Science Shared Resource), NSF grant 1463233, and start-up funds from the U-M Rogel Cancer Center and School of Public Health. JK was supported by NIH grants R01DA048993, R01MH105561 and R01GM124061 and NSF grant IIS2123777. AR and SK were supported by CCSG Bioinformatics Shared Resource 5 P30 CA046592, a gift from Agilent technologies, a Research Scholar Grant from the American Cancer Society (RSG-16-005-01), a Precision health Investigator award from U-M Precision Health to AR along with L. Rozek and M. Sartor, and a University of Michigan MIDAS PODS Grant. SK, VB, and AR were partially supported by the NCI Grant R37-CA214955. SK and AR were also partially supported by The University of Michigan (U-M) startup institutional research funds. NO was supported by the cancer biostatistical training grant #5T32CA083654-18.

Conflict of interest

AR serves as a member for Voxel analytics LLC. and consults for Genophyll, LLC.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2023.1175603/full#supplementary-material

References

Achkar, T., and Tarhini, A. A. (2017). The use of immunotherapy in the treatment of melanoma. J. Hematol. Oncol. 10 (1), 88–89. doi:10.1186/s13045-017-0458-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Agresta, L., KasperHoebe, H. N., and Janssen, E. M. (2018). The emerging role of cd244 signaling in immune cells of the tumor microenvironment. Front. Immunol. 9, 2809. ISSN 1664-3224. doi:10.3389/fimmu.2018.02809

PubMed Abstract | CrossRef Full Text | Google Scholar

Akbani, R., Akdemir, K. C., Aksoy, B. A., Albert, M., Ally, A., Amin, S. B., et al. (2015). Genomic classification of cutaneous melanoma. Cell 161 (7), 1681–1696. doi:10.1016/j.cell.2015.05.044

PubMed Abstract | CrossRef Full Text | Google Scholar

Amgad, M., Sarkar, A., Srinivas, C., Redman, R., Ratra, S., Bechert, C. J., et al. (2019). Joint region and nucleus segmentation for characterization of tumor infiltrating lymphocytes in breast cancer. Med. Imaging 2019 Digit. Pathol. 10956, 129–136.

CrossRef Full Text | Google Scholar

Antohe, M., Nedelcu, R., Nichita, L., Popp, C., Cioplea, M., and Brinzea, A. (2019). Tumor infiltrating lymphocytes: The regulator of melanoma evolution. Oncol. Lett. 17, 1792–1082. ISSN 4155-4161. doi:10.3892/ol.2019.9940

PubMed Abstract | CrossRef Full Text | Google Scholar

Baddeley, A., and Turner, R. (2000). Practical maximum pseudolikelihood for spatial point patterns: (with discussion). Aust. N. Z. J. Statistics 42 (3), 283–322. doi:10.1111/1467-842x.00128

CrossRef Full Text | Google Scholar

Baddeley, A., and Turner, R. (2005). Spatstat: An r package for analyzing spatial point patterns. J. Stat. Softw. 12, 1–42. doi:10.18637/jss.v012.i06

CrossRef Full Text | Google Scholar

Bankhead, P., Fernández, J. A., McArt, D. G., Boyle, D. P., Li, G., Loughrey, M. B., et al. (2018). Integrated tumor identification and automated scoring minimizes pathologist involvement and provides new insights to key biomarkers in breast cancer. Lab. Investig. 98 (1), 15–26. doi:10.1038/labinvest.2017.131

PubMed Abstract | CrossRef Full Text | Google Scholar

Bankhead, P., Loughrey, M. B., Fernández, J. A., Dombrowski, Y., McArt, D. G., Dunne, P. D., et al. (2017). Qupath: Open source software for digital pathology image analysis. Sci. Rep. 7 (1), 1–13.

PubMed Abstract | CrossRef Full Text | Google Scholar

Becht, E., Giraldo, N. A., Lacroix, L., Buttard, B., Elarouci, N., and Petitprez, F. (2016). Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression. Genome Biol. 17 (1), 218. ISSN 1474-760X. doi:10.1186/s13059-016-1070-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Berger, M. F., Hodis, E., Heffernan, T. P., Lawrence, M. S., Protopopov, A., Ivanova, E., et al. (2012). Melanoma genome sequencing reveals frequent prex2 mutations. Nature 485 (7399), 502–506. doi:10.1038/nature11071

PubMed Abstract | CrossRef Full Text | Google Scholar

Bhattacharya, S., Dunn, P., Thomas, C. G., Smith, B., Schaefer, H., Chen, J., et al. (2018). Immport, toward repurposing of open access immunological assay data for translational and clinical research. Sci. data 5 (1), 180015–180019. doi:10.1038/sdata.2018.15

PubMed Abstract | CrossRef Full Text | Google Scholar

Bian, Y., Liu, Y. F., Jiang, H., Meng, Y., Liu, F., Cao, K., et al. (2021). Machine learning for mri radiomics: A study predicting tumor-infiltrating lymphocytes in patients with pancreatic ductal adenocarcinoma. Abdom. Radiol. 46 (10), 4800–4816. doi:10.1007/s00261-021-03159-9

CrossRef Full Text | Google Scholar

Breiman, L. (2001). Random forests machine learning. Mach. Learn. 45, 5–32. View Article PubMed/NCBI Google Scholar. doi:10.1023/a:1010933404324

CrossRef Full Text | Google Scholar

Denkert, C., Von Minckwitz, G., Brase, J. C., Bruno, V. S., Gade, S., Kronenwett, R., et al. (2015). Tumor-infiltrating lymphocytes and response to neoadjuvant chemotherapy with or without carboplatin in human epidermal growth factor receptor 2-positive and triple-negative primary breast cancers. J. Clin. Oncol. 33 (9), 983–991. doi:10.1200/JCO.2014.58.1967

PubMed Abstract | CrossRef Full Text | Google Scholar

Dieci, M. V., Radosevic-Robin, N., Fineberg, S., Van den Eynden, G., Ternes, N., Penault-Llorca, F., et al. (2018). Update on tumor-infiltrating lymphocytes (tils) in breast cancer, including recommendations to assess tils in residual disease after neoadjuvant therapy and in carcinoma in situ: A report of the international immuno-oncology biomarker working group on breast cancer. Seminars cancer Biol. 52, 16–25. doi:10.1016/j.semcancer.2017.10.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Ertosun, M. G., and Rubin, D. L. (2015). Automated grading of gliomas using deep learning in digital pathology images: A modular approach with ensemble of convolutional neural networks. AMIA Annu. Symp. Proc. AMIA Symp. 2015, 1899–1908. ISSN 1942-597X.,

PubMed Abstract | Google Scholar

Franklin, E. L., Roesch, A., Schilling, B., and Schadendorf, D. (2017). Immunotherapy in melanoma: Recent advances and future directions. Eur. J. Surg. Oncol. (EJSO) 43 (3), 604–611. doi:10.1016/j.ejso.2016.07.145

PubMed Abstract | CrossRef Full Text | Google Scholar

Fu, Q., Chen, N., Ge, C., Li, R., Li, Z., Zeng, B., et al. (2019). Prognostic value of tumor-infiltrating lymphocytes in melanoma: A systematic review and meta-analysis. Oncoimmunology 8 (7), e1593806. doi:10.1080/2162402X.2019.1593806

PubMed Abstract | CrossRef Full Text | Google Scholar

Fuchs, N., Meta, M., Schuppan, D., Lutz, N., and Schirmeister, T. (2020). Novel opportunities for cathepsin s inhibitors in cancer immunotherapy by nanocarrier-mediated delivery. Cells 9 (9), 2021. ISSN 2073-4409. doi:10.3390/cells9092021

PubMed Abstract | CrossRef Full Text | Google Scholar

Heindl, A., Nawaz, S., and Yuan, Y. (2015). Mapping spatial heterogeneity in the tumor microenvironment: A new era for digital pathology. Lab. Investig. 95 (4), 377–384. doi:10.1038/labinvest.2014.155

PubMed Abstract | CrossRef Full Text | Google Scholar

Högmander, H., and Aila, S. (1999). Multitype spatial point patterns with hierarchical interactions. Biometrics 55 (4), 1051–1058. doi:10.1111/j.0006-341x.1999.01051.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Idos, G. E., Kwok, J., Bonthala, N., Lynn, K., Gruber, S. B., Qu, C., et al. (2020). The prognostic implications of tumor infiltrating lymphocytes in colorectal cancer: A systematic review and meta-analysis. Sci. Rep. 10 (1), 3360. doi:10.1038/s41598-020-60255-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Kang, J., Johnson, T. D., Nichols, T. E., and Tor, D. W. (2011). Meta analysis of functional neuroimaging data via bayesian spatial point processes. J. Am. Stat. Assoc. 106 (493), 124–134. doi:10.1198/jasa.2011.ap09735

PubMed Abstract | CrossRef Full Text | Google Scholar

Khong, H. T., and Restifo, N. P. (2002). Natural selection of tumor variants in the generation of “tumor escape” phenotypes. Nat. Immunol. 3 (11), 999–1005. doi:10.1038/ni1102-999

PubMed Abstract | CrossRef Full Text | Google Scholar

Khosravi, P., Kazemi, E., Imielinski, M., Elemento, O., and Hajirasouliha, I. (2018). Deep convolutional neural networks enable discrimination of heterogeneous digital pathology images. EBioMedicine 27, 317–328. ISSN 23523964. doi:10.1016/j.ebiom.2017.12.026

PubMed Abstract | CrossRef Full Text | Google Scholar

King, R., Janine, B. I., King, S. E., Nightingale, G. F., and Hendrichsen, D. K. (2012). A bayesian approach to fitting gibbs processes with temporal random effects. J. Agric. Biol. Environ. statistics 17 (4), 601–622. doi:10.1007/s13253-012-0111-0

CrossRef Full Text | Google Scholar

Komura, D., and Ishikawa, S. (2018). Machine learning methods for histopathological image analysis. Comput. Struct. Biotechnol. J. 16, 34–42. doi:10.1016/j.csbj.2018.01.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Q., Wang, X., Liang, F., Yi, F., Xie, Y., Gazdar, A., et al. (2019a). A bayesian hidden potts mixture model for analyzing lung cancer pathology images. Biostatistics 20 (4), 565–581. doi:10.1093/biostatistics/kxy019

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Q., Wang, X., Liang, F., and Xiao, G. (2019b). A bayesian mark interaction model for analysis of tumor pathology images. Ann. Appl. statistics 13 (3), 1708–1732. doi:10.1214/19-aoas1254

CrossRef Full Text | Google Scholar

Lu, Z., Xu, S., Shao, W., Wu, Y., Zhang, J., Han, Z., et al. (2020). Deep-learning–based characterization of tumor-infiltrating lymphocytes in breast cancers from histopathology images and multiomics data. JCO Clin. cancer Inf. 4, 480–490. doi:10.1200/CCI.19.00126

CrossRef Full Text | Google Scholar

Mukherji, B. (2013). Immunology of melanoma. Clin. Dermatology 31 (2), 156–165. ISSN 0738081X. doi:10.1016/j.clindermatol.2012.08.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Negahbani, F., Sabzi, R., Firouzabadi, D., Movahedi, F., Kohandel Shirazi, M., Majidi, S., et al. (2021). Pathonet introduced as a deep neural network backend for evaluation of ki-67 and tumor-infiltrating lymphocytes in breast cancer. Sci. Rep. 11 (1), 1–13.

PubMed Abstract | CrossRef Full Text | Google Scholar

Nguyen, L. T., Pei, H. Y., Nie, J., Liadis, N., Ghazarian, D., Al-Habeeb, A., et al. (2010). Expansion and characterization of human melanoma tumor-infiltrating lymphocytes (tils). PLoS ONE 5 (11), e13940. ISSN 1932-6203. doi:10.1371/journal.pone.0013940

PubMed Abstract | CrossRef Full Text | Google Scholar

Oble, D. A., Robert Loewe, , Yu, P., and Martin, C. M. (2009). Focus on tils: Prognostic significance of tumor infiltrating lymphocytes in human melanoma. Cancer Immun. 9, 3.

PubMed Abstract | Google Scholar

Pagès, F., Galon, J., Dieu-Nosjean, M-C., Tartour, E., Sautès-Fridman, C., and Fridman, W-H. (2010). Immune infiltration in human tumors: A prognostic factor that should not be ignored. Oncogene 29 (8), 1093–1102. ISSN 0950-9232, 1476-5594. doi:10.1038/onc.2009.416

PubMed Abstract | CrossRef Full Text | Google Scholar

Pitcovski, J., Shahar, E., Aizenshtein, E., and Gorodetsky, R. (2017). Melanoma antigens and related immunological markers. Crit. Rev. Oncology/Hematology 115, 36–49. ISSN 10408428. doi:10.1016/j.critrevonc.2017.05.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Ray, M., Kang, J., and Zhang, H. (2015). Identifying activation centers with spatial cox point processes using fmri data. IEEE/ACM Trans. Comput. Biol. Bioinforma. 13 (6), 1130–1141. doi:10.1109/TCBB.2015.2510007

PubMed Abstract | CrossRef Full Text | Google Scholar

Ruifrok, A. C., Johnston, D. A., et al. (2001). Quantification of histochemical staining by color deconvolution. Anal. quantitative Cytol. histology 23 (4), 291–299.

PubMed Abstract | Google Scholar

Saltz, J., Gupta, R., Hou, L., Kurc, T., Singh, P., Nguyen, V., et al. (2018). Spatial organization and molecular correlation of tumor-infiltrating lymphocytes using deep learning on pathology images. Cell Rep. 23 (1), 181–193.e7. doi:10.1016/j.celrep.2018.03.086

PubMed Abstract | CrossRef Full Text | Google Scholar

Seal, S., and Ghosh, D. (2022). Miami: Mutual information-based analysis of multiplex imaging data. Bioinformatics, 38 (15):3818–3826. ISSN 1367-4803, 1460-2059. doi:10.1093/bioinformatics/btac414

PubMed Abstract | CrossRef Full Text | Google Scholar

Seal, S., Vu, T., Ghosh, T., Wrobel, J., and Ghosh, D. (2022). Denvar: Density-based variation analysis of multiplex imaging data. Bioinforma. Adv. 2 (1), vbac039. ISSN 2635-0041. doi:10.1093/bioadv/vbac039

CrossRef Full Text | Google Scholar

Strauss, D. J. (1975). A model for clustering. Biometrika 62 (2), 467–475. doi:10.1093/biomet/62.2.467

CrossRef Full Text | Google Scholar

Sturm, G., Finotello, F., Petitprez, F., Zhang, J. D., Jan, B., Fridman, W. H., et al. (2019). Comprehensive evaluation of transcriptome-based cell-type quantification methods for immuno-oncology. Bioinformatics, 35(14):i436–i445. ISSN 1367-4803, 1460-2059. doi:10.1093/bioinformatics/btz363

PubMed Abstract | CrossRef Full Text | Google Scholar

Van Lieshout, M. N. M., and Baddeley, A. J. (1999). Indices of dependence between types in multivariate point patterns. Scand. J. Statistics, 26(4):511–532. ISSN 0303-6898, 1467-9469. doi:10.1111/1467-9469.00165

CrossRef Full Text | Google Scholar

Wang, S., Wang, T., Yang, L., Yang, D. M., Fujimoto, J., Yi, F., et al. (2019). Convpath: A software tool for lung adenocarcinoma digital pathological image analysis aided by a convolutional neural network. EBioMedicine 50, 103–110. doi:10.1016/j.ebiom.2019.10.033

PubMed Abstract | CrossRef Full Text | Google Scholar

Weaver, J. D., Stack, E. C., Buggé, J. A., Hu, C., McGrath, L., Mueller, A., et al. (2022). Differential expression of ccr8 in tumors versus normal tissue allows specific depletion of tumor-infiltrating t regulatory cells by gs-1811, a novel fc-optimized anti-ccr8 antibody. OncoImmunology 11 (1), 2141007. ISSN 2162-402X. doi:10.1080/2162402X.2022.2141007

PubMed Abstract | CrossRef Full Text | Google Scholar

Weiss, S. A., Han, S. W., Lui, K., Tchack, J., Berman, R., et al. (2016). Immunologic heterogeneity of tumor-infiltrating lymphocyte composition in primary melanoma. Hum. Pathol. 57, 116–125. ISSN 00468177. doi:10.1016/j.humpath.2016.07.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu-Sung, S., Masanao, Y., and Yu-Sung, M. (2015). Package ‘r2jags’. R package version 0.03-08. Available at: http://CRAN.R-project.org/package=R2jags.

Google Scholar

Yuan, Y., Failmezger, H., Oscar, M. R., Ali, H. R., Gräf, S., Chin, S-F., et al. (2012). Quantitative image analysis of cellular heterogeneity in breast tumors complements genomic profiling. Sci. Transl. Med. 4 (157), 157ra143. doi:10.1126/scitranslmed.3004330

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y., and Zhang, Z. (2020). The history and advances in cancer immunotherapy: Understanding the characteristics of tumor-infiltrating immune cells and their therapeutic implications. Cell. Mol. Immunol., 17(8):807–821. ISSN 1672-7681. doi:10.1038/s41423-020-0488-6

PubMed Abstract | CrossRef Full Text | Google Scholar