Ethics statement

The study was approved by the Yale University Human Investigation Committee protocols #9505008219 and #1608018220; or local institutional protocols which approved the patient consent forms or, in some cases waiver of consent when retrospectively collected archive tissue was used in a de-identified manner.

Overview

We start by presenting the problem setting. Consider n data points $$ {\left\{x_i\right\}}_{i=1}^n$$ from a hidden manifold $$ \mathcal{M}$$ embedded in a high-dimensional Euclidean space $$ {\mathbb{R}}^k$$. Assume we do not have direct access to these data points; instead, these data points are measured through m observation functions $$ f_j:\mathcal{M}\to {\mathbb{R}}^d$$, where j = 1, …, m is the index of the observation function. Given n multi-feature observations f_j(x_i) of the data point x_i for i = 1, …, n, each consisting of m features j = 1, …, m, our goal is to recover the data points x_i on the hidden manifold $$ \mathcal{M}$$.

In the context of IMC, the data points represent the treatment outcome based on n spatial positions located at ROIs within pathology slides of tissues from several patients. At each spatial position i = 1, …, n, the observations f_j(x_i) for j = 1, …, m are the expression levels of m markers. Each observation $$ f_j\left(x_i\right)\in {\mathbb{R}}^d$$ is a patch of d pixels of the expression level image of marker j at position i.

To simplify the presentation of our approach, we begin with an illustrative localization problem, which is simpler than the IMC problem. Suppose we have a surface $$ \mathcal{M}$$ and objects located at x_i on $$ \mathcal{M}$$. The locations x_i are hidden, but measured through m sensors, such that for each location x_i we have m multi-feature observations $$ {\left\{f_j\left(x_i\right)\right\}}_{j=1}^m$$. That is f_j(x_i) is a d-dimensional observation of sensor j when the object is at x_i. The goal is to recover the locations x_i on the surface $$ \mathcal{M}$$ given $$ {\left\{f_j\left(x_i\right)\right\}}_{j=1}^m$$.

Our approach consists of two stages. At the first stage, we construct a graph for each data point x_i in order to capture associations between its m multi-feature observations $$ {\left\{f_j\left(x_i\right)\right\}}_{j=1}^m$$. Capturing such mutual-relationships is natural in the context of localization problems, e.g., the triangulation property [18] in which the relative locations of the sensors are exploited. In addition, these mutual-relations are typically more robust to noise in comparison with the nominal values of the multi-feature observations, f_j(x_i), themselves. Concretely, consider the m observations $$ {\left\{f_j\left(x_i\right)\right\}}_{j=1}^m$$ associated with the data point x_i. Each observation f_j(x_i) forms a single node in the graph, hereby denoted as node j, giving rise to a graph with a total of m nodes. The graph we consider is the complete graph, where the weights of the edges are determined based on the Euclidean distance between the corresponding observations: the weight of the edge connecting nodes j and k is proportional to exp{−‖f_j(x_i) − f_k(x_i)‖²}. Then, we compute the SSD of a random walk defined on this graph at each location. SSD has a vector form that embodies the multi-feature inter-relationships of the data point x_i.

At the next step, we define a second graph based on the SSDs, characterizing the points $$ {\left\{x_i\right\}}_{i=1}^n$$. Concretely, each data point x_i is represented by a node, and the pairwise distances between the SSDs form the adjacency matrix of the graph. Then, we apply a particular manifold learning technique, diffusion maps [4], to this graph. This application facilitates the recovery of the hidden manifold $$ \mathcal{M}$$ in the sense that an embedding of the points x_i is learned, such that the distances between the embedded points respect a notion of an intrinsic distance (the diffusion distance [4]) on $$ \mathcal{M}$$. The application of diffusion maps to the second graph in the context of the localization problem gives rise to an embedding that serves as an accurate representation of the hidden locations of the data point. In Localization toy problem, we demonstrate the proposed method on several simulations of localization toy problems.

We remark that the IMC problem and the localization problem share many aspects. For example, in both problems, the multi-feature observations are noisy and the mutual-relationship between them carry important information. Yet, there is a particular aspect that makes the IMC problem more challenging; while the points on the hidden manifold in a localization problem are homogeneous because they all represent location coordinates, the points in the IMC problem could be significantly different due to the large variability in the tissue structure. Importantly, the proposed method accommodates the joint processing of such different points through their representation by the SSD.

Proposed method

The first step of the proposed method is to construct an undirected weighted graph $$ {\mathcal{G}}_i=({\mathcal{V}}_i,{\mathcal{E}}_i,{\mathit{W}}_i)$$ for each data point $$ x_i\in \mathcal{M}\subset {\mathbb{R}}^k$$, where the vertex set is $$ {\mathcal{V}}_i=\{f_1\left(x_i\right),f_2\left(x_i\right),...,f_m\left(x_i\right)\}$$, the edge set is $$ {\mathcal{E}}_i\subseteq {\mathcal{V}}_i\times {\mathcal{V}}_i$$, and the graph weights matrix $$ {\mathit{W}}_i\in {\mathbb{R}}^{m\times m}$$ is given by the Gaussian kernel

Next, we define a random walk on the graph $$ {\mathcal{G}}_i$$. Let $$ {\mathit{P}}_i\in {\mathbb{R}}^{m\times m}$$ be a row stochastic matrix given by

The transition probability matrix P_i is self-adjoint and compact, and therefore, the spectral decomposition of $$ {\mathit{P}}_i^t$$ for t > 0 is given by

Interestingly, in the special case where t = 1, the probability of the node f_k(x_i) to stay in place is given by

Note that $$ {\mathit{\varphi }}_1\in {\mathbb{R}}^m$$ is the left eigenvector of P_i corresponding to λ₁ = 1, satisfying:

We will use π_i as a new characteristic vector, or a signature, of x_i, and consequently, the induced pairwise distances ||π_i − π_i′||, where i, i′ ∈ 1, …, n, will be used as the desired distances between the respective graphs $$ {\mathcal{G}}_i$$ and $$ {\mathcal{G}}_{i^{\prime }}$$ for recovering $$ \mathcal{M}$$. At first glance, using π_i may seem too simplistic. Instead, one could use the broad spectral information. Consequently, define the Diffusion Kernel Signature (DKS) by

Since λ_j is in descending order and in [0, 1], the weight they assign to the eigenvectors in (12) becomes smaller as t increases. As a result, the DKS can be viewed as a low-pass filter, which controls the spectral bandwidth. In addition, the DKS can be recast in terms of the diffusion distance, a notion of distance induced by diffusion maps [4] that was shown useful in a broad range of applications, e.g., [20–22]. For more details see Diffusion maps. Specifically, when t = 1, we can show that

We note that DKS has already appeared in previous work in the context of spectral distances in [23, 24], where it was shown that it describes the underlying geometry of $$ \mathcal{M}$$. We show in the following that the seemingly simple SSD, despite the lack of broad spectral information as in the DKS, still carries substantial information.

Note that ϵ is a scale parameter of the Gaussian kernel, where it can be used to infer locality. If ϵ is set to a small value, then π_i captures local properties. Conversely, if ϵ is large, then π_i represents the global structure. As a result, a multiscale signature can be formed, consisting of multiple SSDs π_i computed with different values of ϵ.

The final stage of our method is building a low-dimensional representation of all the data points $$ {\left\{x_i\right\}}_{i=1}^n$$. To this end, we apply diffusion maps to the corresponding characteristic vectors (signatures) $$ {\left\{{\mathit{\pi }}_i\right\}}_{i=1}^n$$ as follows. First, we build a global graph G⁽²⁾ whose nodes are π_i and edge weights are determined by a Gaussian kernel based on the l₁ distance between π_i.

That is, the global graph weights matrix W⁽²⁾ is defined by

Second, we construct a random walk, denoting its transition probability matrix by P⁽²⁾. Third, we apply the eigendecomposition to P⁽²⁾. Fourth, we set the dimension of the new representation according to the variant of the Jackstraw method [6], as shown in Determining the dimension of data.

The entire method is summarized in Box 1 and a block diagram is illustrated in Fig 1.

Fig 1

Illustrative diagram of the proposed method.

(a) For each data point x_i, we build a local graph $$ {\mathcal{G}}_i$$ based on its multi-feature observations $$ {\left\{f_j\left(x_i\right)\right\}}_{j=1}^m$$. (b) We construct a random walk with transition probabilities matrix P_i on $$ {\mathcal{G}}_i$$. (c) We extract the SSD signature π_i from P_i. (e) We collect the SSDs of $$ {\left\{x_i\right\}}_{i=1}^n$$ into an SSD representation matrix. (g) Subsequently, the matrix is subjected to a nonlinear dimensionality reduction using diffusion maps by the construction of the global graph G⁽²⁾ and the corresponding random walk with P⁽²⁾. (f) Via eigenvalue decomposition, we obtain a low-dimensional representation Ψ(x_i) for i = 1, …, n, which is used in the subsequent tasks (g).

Theoretical analysis

We propose a statistical model that allows for a tractable analysis, showing the advantages of the SSD signature. Consider a data point $$ x_i\in \mathcal{M}$$ and denote the set of the observations by $$ S={\left\{f_j\left(x_i\right)\right\}}_{j=1}^m$$. Assume the j-th observation $$ f_j\left(x_i\right)\in {\mathbb{R}}^d$$ is a realization of a d-dimensional random vector $$ V_j^i$$ following a multivariate normal distribution given by

By collecting the m observations $$ {\left\{f_j\left(x_i\right)\right\}}_{j=1}^m$$ and denoting the covariance matrix between the k and l observation functions by $$ {\mathbf{\Sigma }}_{k,l}^i$$, we obtain an md-dimensional vector that can be viewed as a realization of the random vector $$ V^i=(V_1^i,\dots ,V_m^i)$$ with the corresponding multivariate Gaussian distribution given by

Definition 1 (empirical mean) Given a set Γ and some real function on the set $$ \mathit{q}\in {\mathbb{R}}^{|\Gamma |}$$, and a subset Ω ⊂ Γ, the empirical mean of q in Ω is defined by

Definition 2 (heterogeneity) Define the heterogeneity of a data point x_i by

The heterogeneity $$ {\mathit{h}}_i\in {\mathbb{R}}^m$$ captures the mutual-relationships between the expected values of the observations $$ {\mu }_j^i$$ of a particular data point x_i; if Θ(j)ⁱ significantly deviates from μⁱ, then h_i(j) is large. Conversely, if Θⁱ(j) is close to μⁱ, then h_i(j) is close to zero.

Definition 3 (weighted heterogeneity) Define g_i as a weighted heterogeneity, whose jth element is given by

Under the considered statistical model, with the above definitions, the SSD π_i can be written explicitly.

Proposition 1 The j-th element of π_i can be approximated by

The derivation is based on the Taylor expansion of the Gaussian function in Eq (1), where ϵ is the scale of the function. The proof appears in S1 Appendix.

In order to give some intuition, we consider the following special cases, where the SSD assumes a simpler form.

Special case 1

Suppose that the random vectors of the observation functions are independent and identically distributed, i.e., $$ {\mathbf{\Sigma }}_{k,l}^i=\mathbf{0}$$ ∀l ∈ {1, …, m} and k ≠ l. In this case, the k-th element of π_i is

where

Note that a small value is assigned to π_i(k) if the weighted heterogeneity g_i(k) is large. In contrast, a large value is assigned to π_i(k) if g_i(k) is small. As a consequence, π_i(k) carries the heterogeneity information of the observations.

Special case 2

When the kernel scale ϵ → ∞, the information about the heterogeneity of each observation is lost, since the same weights are assigned to all the edges. As a consequence, π_i becomes just a constant vector, given by

where $$ \mathbf{1}\in {\mathbb{R}}^m$$ is an all-ones vector.

Binary hypothesis testing

Suppose that $$ {\left\{x_i\right\}}_{i=1}^n\in \mathcal{M}$$ are realizations of a random variable X, which follows a bimodal distribution stemming from two hypotheses: $$ {\mathcal{H}}_1$$ and $$ {\mathcal{H}}_2$$; $$ {\mathcal{H}}_1$$ has probability α and $$ {\mathcal{H}}_2$$ has probability (1 − α), where 0 < α < 1. Denote the set of data points from hypothesis $$ {\mathcal{H}}_1$$ by Ω₁ and the set of data points from hypothesis $$ {\mathcal{H}}_2$$ by Ω₂. Recall that for each data point x_i, f_j(x_i) are the realizations of the elements of the random vector Vⁱ. Since Vⁱ depends on the random variable X, assume that Vⁱ also follows a bimodal distribution, which is induced by the bimodal distribution of X. Particularly, consider a Gaussian setting, where Vⁱ is sampled from $$ \mathcal{N}({\mathit{m}}_1,{\mathit{\Sigma }}_1)$$ with probability α and from $$ \mathcal{N}({\mathit{m}}_2,{\mathit{\Sigma }}_2)$$ with probability of (1 − α). Similarly, assume that the random variable $$ V_j^i$$ follows a bimodal distribution: sampled from $$ \mathcal{N}({\mu }_j^1{\mathbf{1}}_d,{\sigma }_j^1{\mathit{I}}_d)$$ with probability α and from $$ \mathcal{N}({\mu }_j^2{\mathbf{1}}_d,{\sigma }_j^2{\mathit{I}}_d)$$ with probability (1 − α), where the respective probability density functions are $$ f\left(v\right|{\mu }_j^1,{\sigma }_j^1)$$ and $$ f\left(v\right|{\mu }_j^2,{\sigma }_j^2)$$.

A naïve approach for binary hypothesis testing would be to directly compare the densities of the two hypotheses for each observation separately. Particularly, based on the realizations from only one observation function f_j, the average probability of error attained in a Bayesian setting with the MAP estimator [5] is given by

According to [25], consider $$ \frac{|{\sigma }_j^1-{\sigma }_j^2|}{\sigma \left(j\right)}\le \frac{2}{3}$$, where $$ \sigma \left(j\right)=max\{{\sigma }_j^1,{\sigma }_j^2\}$$. The total variation above between two Gaussian distributions is bounded by

We seek another more discriminative approach for binary hypothesis testing. For this purpose, we propose a method based on the SSDs. Since the obtained SSDs represent probability distributions, the average probability of error is given by

According to Proposition 1, the total variation between two SSDs associated with data points from two hypotheses can be explicitly expressed as the l₁ distance given by

The total variation of the measurements in Eq (29) and the total variation between the SSDs in Eq (32) can be used to distinguish between the two hypotheses. In the following, we specify the conditions, under which the total variation based on the SSDs in Eq (31) is larger, and hence, leading to smaller error compared to the standard MAP estimator using a single observation as specified in Eq (28).

Proposition 2 Suppose $$ {\mu }_j^1={\mu }_j^2$$ and $$ {\sigma }_j^1={\sigma }_j^2$$, which imply by definition (or by Eq (30)) that the total variation between the distributions corresponding to the two hypotheses is zero, i.e.,

This means that not only the standard MAP estimator but also any estimator based directly on single channel observations cannot distinguish between the two hypotheses. Conversely, from Eq (32), the SSDs may carry a distinction capability, that is,

Proposition 2 demonstrates that there are cases where a single observation cannot be used for distinguishing between the two hypotheses. However, in such cases, the SSDs may enable us to distinguish the hypotheses due to possible differences in either the heterogeneity or the covariances. Proposition 2 is further demonstrated in the context of the localization toy problem in Simulation 2.

To further expand the analysis, we make the following assumptions.

Assumption 1 The empirical mean of the weighted heterogeneity is approximately the same under the two hypotheses:

Assumption 2 The empirical mean of the covariance matrices is approximately the same under the two hypotheses:

Note that if Assumptions (A.1) and (A.2) hold, implying that $$ 2ϵm-2m({\langle {\mathit{g}}^2\rangle }_{{\Omega }_1}+{\langle \mathbf{\Sigma }\rangle }_{{\Omega }_1})\simeq 2ϵm-2m({\langle {\mathit{g}}^2\rangle }_{{\Omega }_2}+{\langle \mathbf{\Sigma }\rangle }_{{\Omega }_2})$$, then the l₁ distance between the SSDs in Eq (32) can be recast as

Proposition 3 Suppose that Assumptions (A.1) and (A.2) hold. Suppose $$ {\mu }_j^1={\mu }_j^2$$, which implies that the upper bound of the total variation at the j-th element in Eq (30) only depends on the variance

In addition, suppose that (i) the covariance between the jth observation and the other observations under the two hypotheses is approximately equal $$ {\langle {\mathbf{\Sigma }}_{k,·}\rangle }_{{\Omega }_1}\simeq {\langle {\mathbf{\Sigma }}_{k,·}\rangle }_{{\Omega }_2}$$, (ii) the empirical mean of the difference of variance and weighted heterogeneity of the two hypotheses is greater than the difference of j-th variance, i.e., $$ \langle {({\mathit{\sigma }}^1-{\mathit{\sigma }}^2)}^{\top }({\mathit{g}}_1^2-{\mathit{g}}_2^2)\rangle \ge |{\sigma }_j^1-{\sigma }_j^2|$$, and (iii) the weighted heterogeneity of $$ {\mathcal{H}}_1$$ is sufficiently large such that $$ {\langle {\mathit{g}}^2\rangle }_{{\Omega }_1}\ge ϵ-{\langle \mathbf{\Sigma }\rangle }_{{\Omega }_1}-\sigma \left(j\right)$$. Then, the l₁ distance between the SSDs can be recast and bounded from below by

It follows that

This proposition implies that when the assumptions hold, the probability of error based on SSD, which indirectly takes into account the mutual-relations between all observations, facilitates a better distinction of the two hypotheses compared to the standard MAP estimator computed from the best sensor. This property is further demonstrated in the localization toy problem in Simulation 3.

Proposition 4 Suppose the conditions of Proposition 3 hold. In addition, suppose that the random vectors of the observations are independent and identically distributed, i.e., $$ {\sigma }_{k,j}^1={\sigma }_{k,j}^2=0$$ ∀k ≠ l, then

This proposition shows that the SSDs enable us a better distinction between the two hypotheses compared to the MAP estimator based on the distributions of the multi-feature observations. In other words, the heterogeneity comprising the SSD has a significant contribution to the ability to recover the information about the latent data points $$ {\left\{x_i\right\}}_{i=1}^n$$ on the underlying manifold $$ \mathcal{M}$$, thereby leading to accurate binary hypothesis testing.

Imaging mass cytometry (IMC)

IMC is a relatively new imaging method, which enables to examine tumors and tissues at subcellular resolutions, giving rise to images consisting of the intensities of multiple proteins [13–15]. This acquisition platform, combined with computational methods, has recently been the subject of many studies. Various image processing and analysis techniques for IMC datasets can be found in [26], where it is shown that single-cell segmentation can be accomplished successfully with supervised classifiers, leading to the characterization of cell co-occurrence and cell composition of different types of tissues and samples. In [27], an IMC dataset with 37 markers is used for cell segmentation and cell clustering based on random 125 × 125μm² patches collected from breast cancer patients. This dataset is jointly analyzed with multi-platform genomics data, where it is shown that classifiers can be iteratively trained in a supervised manner to learn from the IMC pixels the corresponding cell expression levels. In [14], spanning-tree progression analysis combined with samples’ type provided by pathologists is used for cell population and cell transition identification. In contrast to these methods, our approach focuses on extracting the mutual-relationships between markers at likely tumor cells regions at large, circumventing cell segmentation.

In this work, our goal is to identify the sensitivity of lung cancer subjects to treatment with PD-1 axis blockers, given their IMC multiplexed observations. More specifically, we aim at a binary prediction task: identifying whether the subjects responded or did not respond to the treatment. We analyze two IMC datasets consisting of baseline/treatment tumor samples from non-small cell lung cancer subjects profiled with 29 markers, representing phenotype and functional properties of both tumor and immune cells. The markers are denoted by LipoR, VIM, T-BET, CD47, Cytokeratin, CD45RO, PD-L1, GAPDH, B7-H3, LAG-3, TIM-3, FOXP3, CD4, B7-H4, CD68, PD1, CD20, CD8, CD25, VISTA, KI-67, B2M, CD3, IDO-1, PD-L2, GZB, Histone 3, DNA1 and DNA2. The resolution of the IMC images is 1 μm² per pixel. The first dataset, denoted Dataset 1, consists of 55 subjects (samples), and the second dataset, denoted Dataset 2, consists of 29 subjects (samples). These subjects received treatment with PD-1 axis blockers. Based on the clinical sensitivity to the treatment, the subjects are categorized as: durable clinical benefit (denoted with the label responders) and no durable benefit (denoted with the label non-responders).

Prior to our analysis, a standard pre-processing using z-score normalization is applied to each marker. We remark that the mean and the standard deviation are computed based only on pixels in which the marker has non-zero values. In addition, we apply a 3 × 3 median filter to every image in Dataset 2. We note that the median filter is not applied to Dataset 1, because the typical expression levels in this dataset are sparse, and in this situation, application of a filter may destroy the signal. This is demonstrated in S1 Fig, where we apply the median filter to the expression levels of few important markers and observe a degenerate result. The markers exhibiting such sparse expression levels in Dataset 1 are VIM, T-BET, CD45RO, PD-L1, GAPDH, B7-H3, LAG-3, FOXP3, B7-H4, PD1, CD20, CD8, CD25, VISTA, KI-67, CD3, PD-L2 and GZB. Conversely, in Dataset 2, the typical expression levels are dense, and as a result, the median filter enhances the content of the images.

Our analysis does not consider the entire image, but rather focuses on ROIs located at the highest Cytokeratin expression levels, as Cytokeratin is expressed only in the tumor cells. As noted above, this circumvents cell segmentation that is typically required in other analyses techniques. At each ROI, we consider a “stack” of m image patches from all the m markers, where each patch consists of d = b × b pixels. We select N ROIs per sample by searching the patches with the maximal mean value of Cytokeratin. This is implemented by a 2D convolution applied to the Cytokeratin image with a constant kernel of size b × b. We deliberately avoid patch overlap by assigning zero values to the pixels of each ROI, once it is selected.

Using the present work notation, the intrinsic representation of the information embodied at each ROI is denoted by x_i, which is assumed to be a data point from some hidden manifold $$ \mathcal{M}$$. Our working hypothesis is that the distribution of $$ {\left\{x_i\right\}}_{i=1}^n$$ on the hidden manifold $$ \mathcal{M}$$ is bimodal, which is induced by the sensitivity of the subjects to the treatment; data points x_i at ROIs within tissues of responders are located in one region of the manifold, and data points x_i at ROIs within tissues from non-responders are located in another region of the manifold. Given a data point x_i, our two hypotheses, $$ {\mathcal{H}}_{\text{r}}$$ and $$ {\mathcal{H}}_{\text{n}}$$, are whether x_i is a realization from the distribution of responders or non-responders, respectively. Here, n = N × P is the total number of ROIs across all the samples, where P is the number of samples and N is the number of ROIs we consider in each sample. Next, recall that the intrinsic representation x_i is hidden. Instead, the accessible observations are the expression levels of the markers at the ROIs, which are represented mathematically by the observation functions $$ f_j\left(x_i\right)\in {\mathbb{R}}^d$$ for j = 1, …, m, where m = 29 is the number of markers. The domain of the observation functions is the hidden intrinsic manifold $$ \mathcal{M}$$ and the range is of dimension d = b × b, which is the size of a patch of an image of one marker. Namely, the values assigned by f_j(x_i) to x_i correspond to the expression level of marker j at the ROI associated with x_i. The set of patches of all the markers at a specific ROI x_i is a set of multi-feature observations $$ {\left\{f_j\left(x_i\right)\right\}}_{j=1}^{29}$$. In Fig 2, we illustrate the IMC multiplexed observations from a single subject and the set up under consideration.

Fig 2

An illustration of the IMC multi-feature observations from a single subject.

In the IMC data of a single subject, we focus first on the Cytokeratin marker, which is used as an indicator of tumor cells. We select N ROIs located at the highest Cytokeratin expression levels. These ROIs are patches of size b × b μm², where here b set to 35. We assume that these ROIs have intrinsic hidden states represented by x_i. The expression of all 29 markers at these ROIs are viewed as the multi-feature observations $$ {\left\{f_j\left(x_i\right)\right\}}_{j=1}^{29}$$ for i = 1, …, N.

The identification of the subject’s response to treatment from the IMC data is based on the application of the proposed method described in Box 1. This algorithm fits well the problem at hand due to the following main reasons. First, directly comparing observations from the different markers, f_j(x_i), is inapplicable since each ROI comprises different cells and different tissue structures. Our method circumvents this problem by computing an intrinsic signature of the observations. Second, due to the different dynamics of the nominal values of the observations from the different markers, a naïve concatenation of the multi-feature observations is inadequate (in contrast to the localization example).

Before applying the proposed method, we test empirically that Assumptions (A.1) and (A.2) hold. Note that for this test only, the true response status is used. To test Assumption (A.1), we compute

We analyze the two datasets separately, since the datasets were collected around a year apart, and internal acquisition system parameters were modified during that time. We apply the proposed method presented in Box 1 to the observations, $$ {\left\{f_j\left(x_i\right)\right\}}_{j=1}^{29}$$ for i = 1, …, n, resulting in a low-dimensional representation of the ROIs. Then, we apply to the low-dimensional representation an RBF SVM classifier with a leave-one-subject-out (LOSO) cross-validation [5] in order to predict the response to treatment. In order to assess the prediction performance for each subject, we compute the average of the prediction results of all the ROIs of that subject. We note that the prediction is based on features computed per patch (rather than per subject). Therefore, in practice, the number of samples used for the cross-validation is 55 × N for Dataset 1 and 29 × N for Dataset 2. Importantly, at each cross-validation iteration, all N patches of a subject were removed from the training set, and were only used for testing the classifier.

We compare the SSD-based representation obtained by the propose method to other representations obtained by three competing algorithms. The first is a direct application of diffusion maps to the sets of multi-feature observations $$ {\left\{f_j\left(x_i\right)\right\}}_{j=1}^{29}$$ for i = 1, …, n. The second and the third are based on HKS [16] and WKS [17], respectively, replacing the SSD as the features of each ROI at the first stage of the proposed method in Box 1. The dimensions of the representation obtained by the proposed method and by the three competing algorithms are determined by a variant of the Jackstraw method, described in Determining the dimension of data.

Fig 3 presents the predictions of treatment response by the SSD, HKS, WKS based algorithms as well as DM. We show the 3D t-SNE visualization [28] of the patch representation obtained as an output of the different algorithms and the confusion matrix of the prediction obtained by an RBF SVM classifier applied to the patch representation. To complement the results, in Table 1, we present the area under the ROC curve (AUC) of the treatment response predictions. We note that the AUC obtained for Dataset 2 without the median filter pre-processing is 0.931. For each algorithm, the presented prediction results are based on the patch size and number of patches configuration that yielded the best empirical performance using cross-validation as described above. The best configuration (patch size b × b and number of patches per subject N) of each algorithm is presented in Fig 3 on the left.

Table 1

ROC AUC for predictions of treatment sensitivity.

	Proposed Method	DM	HKS	WKS
Dataset 1	0.9455	0.8364	0.7455	0.7818
Dataset 2	1	0.8276	0.9655	0.9310

Fig 3

Treatment response predictions.

(A) Prediction for Dataset 1. (B) Prediction for Dataset 2. In each panel, at each row, we plot the best algorithm configuration, the 3D t-SNE visualization of the patches representation colored by the response status, and the confusion matrix of the response prediction obtained by a leave-one-subject-out cross-validation using an RBF SVM classifier. The rows present the results of the different methods.

In the t-SNE plots, we observe that the unsupervised separation of the patches from responders and non-responders is most pronounced in the low-dimensional representation obtained by the proposed method. This distinct visual separation, which was obtained by the proposed algorithm without access to any response outcome information. The color in the figures indicates the response status. Note that the embedding is obtained by the unsupervised algorithm and the color labels are overlaid to demonstrate the degree of separability between patches from responders and non-responders, which implies that this patch representation is informative and useful for the subsequent response prediction. This result, which is obtained in an unsupervised manner, distinguishes the current work from the computational methods for IMC data described above that rely on supervised analysis. Indeed, we observe that the prediction accuracy obtained based on the proposed method is superior compared to the other three competing methods. We note that Dataset 1 consists of 26 responders and 29 non-responders, so that the chance level of accurate prediction of a subject with durable clinical benefit is 47.27%, and Dataset 2 consists of 11 responders and 18 non-responders, thus, the chance level of accurate prediction is 37.93%.

The derivations in Theoretical analysis imply that the capability to distinguish between the two hypotheses (sensitivity or insensitivity to treatment) highly depends on the mutual relationships between the markers, which we explicitly define and term heterogeneity (Definition 2). Since our empirical study demonstrates that our approach facilitates a distinct separation between ROIs of subjects according to their treatment response, by the theoretical analysis, we conclude that the heterogeneity between the markers is where the information about the sensitivity to treatment lies.

We examine the sensitivity of the tested algorithms to the choice of the hyperparameters: the number of ROIs per subject N and the size of the patch b × b. In S2 Fig, we plot heatmaps of the treatment prediction accuracy obtained based on different choices of hyperparameters. We observe that within a relatively wide range of parameter values, the performance of the proposed method in Box 1 is high and insensitive to the particular choice of parameters. In addition, we note that the range of patch size where high performance is attained is centered at size 45 × 45 μm². As the size of tumor cells vary within a range of 10–30 μm² (a mean of 20 μm²), and the size of a lymphocyte ranges from 8–12 μm² (a mean of 10 μm²), it implies that high performance is attained when patches are likely including more than one cell and cell type. Conversely, we observe that the competing methods do not show the same degree of robustness to the choice of hyperparameters, and good performance is attained only for very particular (isolated) parameter values.

To further exploit the robustness of the proposed method, we implemented an ensemble of the classifiers based on different values of hyperparameters. The implementation of the combination is based on [29]. The performance of the combined classifiers is presented in S3 Fig. We observe that the prediction accuracy of this ensemble is comparable to the prediction accuracy obtained based on the classifier with the best parameter configuration, demonstrating that the particular choice of hyperparameters can be circumvented.

Our premise is that the resulting high prediction accuracy is attributed mainly to the unsupervised informative representation obtained by our approach, rather than the classifier type. To support this claim, we repeat the analysis by replacing the RBF SVM classifier with Random Forest [5]. The results are shown in S4 Fig and demonstrate comparable performance and similar trends.

To show that Stage 2 of the proposed method in Box 1 is essential, we plot in S5 Fig heatmaps of the multi-feature observations and the SSD features resulting from Stage 1 of the algorithm. The heatmaps of the multi-feature observations demonstrate that there is no obvious difference between responders and non-responders, implying that the response status prediction is a non-trivial task. In addition, inspection of the heatmaps with the SSD features does not reveal apparent distinction between responders and non-responders. Recalling that after Stage 2, as presented in the t-SNE plots in Fig 3, this distinction becomes evident, demonstrating the contribution of Stage 2.

Discussion

We presented a two-step graph analysis approach. The first step is applied to the multi-feature observations of the data points, where their mutual-relationships are extracted. This step is implemented by computing the SSD of a random walk defined on the graph whose nodes are the observations. The resulting SSD can be viewed as a signature or a characteristic vector of the data point and is analogous to traditional signatures from other domains, such as the heat kernel signature (HKS) [16] or the wave kernel signature (WKS) [17] in the field of shape analysis, and PageRank [30] in web page ranking (see Related work). The second step is applied to the signatures obtained at the first step, for the purpose of constructing an intrinsic low-dimensional representation of all the data points.

Previous attempts to analyze such mulitplexed datasets involve various approaches, including direct comparisons of the marker expressions, the cell morphology, and interactions in cell neighborhoods, to name but a few [26]. Our method introduces a new approach, building a new representation of the multiplexed data in two steps. Since each of the steps involves a construction of a graph, the entire procedure can be viewed as building a graph of graphs.

While the algorithm is described in a general setting of multiplexed data fusion, our theoretical analysis is focused on binary hypothesis testing. In comparison with a traditional statistical estimation approach, we show that our method exhibits advantages, implying that the mutual-relationships between the multi-feature observations are well captured. In the context of IMC, this could minimize the effect of deviation in individual marker scores and cell/tissue heterogenenity.

We apply the proposed method to two IMC datasets and show that solely from the imaging data, we can distinguish between two different sensitivity levels to treatment. Since our approach does not rely on rigid prior knowledge or access to labels, it has the potential of identifying biological relevance of novel parameters or marker patterns in treatment responses by analyzing dominant factors contributing more to the model stratification. Importantly, we remark that in contrast to common practice, the proposed approach does not require cell-segmentation as a precursor.

In addition to the demonstrated advantage of the proposed method over the competing methods in terms of superior prediction accuracy of treatment response, the proposed method is also more computationally efficient. In S6 Fig, we present the run time of the two stages of the proposed method in Box 1 and compare them to the run time of the three competing methods. We note that the HKS and WKS replace the proposed SSD in Stage 1, but, their respective Stage 2 are similar. In Stage 1, we observe that the proposed method in Box 1 is faster than the competing methods. The main difference between the algorithms in Stage 1 is due to the fact that the SSD is proportional to the degree of the local graph, and as a result, it can be computed efficiently from the degree vector. Conversely, both HKS and WKS require eigenvalue decomposition, which is computationally more demanding. In addition, the direct application of DM just involves Stage 2 of the proposed algorithm. Seemingly, the run time of the algorithms in Stage 2 should have been the same, yet, we observe that DM is slower. This difference is attributed to the significantly different size of the feature space. In DM, the multi-feature observations are simply concatenated, giving rise to feature space (input of Stage 2) of size (b × b × m) × (P × N), where b × b is the size of patch, m is the number of biomarkers, and P is the number of subjects. Conversely, in Stage 2 of Algorithm 1, HKS, and WKS, the feature space (input of Stage 2) is only of size m × (P × N).

It is conceivable that the most important hyperparameter of our method is the scale parameter. In S7 Fig, we present a toy example demonstrating that different values of the scale parameter ϵ lead to multiscale signatures capturing local and global features. We demonstrate that different scales facilitate the extraction of different features of the data. In future work, we plan to further explore the role of the scale and to devise multiscale signatures. Another possible direction for future research relies on the fact that our method is general and can be extended to other multiplexed datasets. For example, hyper spectral imaging, sensor networks, spatial multiplexed proteomics, and spatial transcriptomics assays is a representative subset of distinct technologies from diverse domains of science and engineering that share common data structures. The data in all these modalities consist of high-dimensional multivariate observations (m-dimensional feature space) collected at different spatial positions, and therefore, can be analyzed using similar computational methodologies. Furthermore, in many studies practitioners collect datasets consisting of multiple spatial assays of this type, each capturing such data from a single biological sample, patient, or hyper spectral image, etc. Each of these spatial assays could be characterized by several regions of interest (ROIs), giving rise to a setting similar to the IMC problem considered here. Specifically, we plan to examine applications of the proposed graph of graphs analysis to spatial transcriptomics such as Slide-seq [31], High-Density Spatial Transcriptomics [32, 33], MIBI-TOF [34], and DBiT-seq [35].

Diffusion maps

Manifold learning is a class of nonlinear techniques that embeds high dimensional data points into a low dimensional space, relying on the assumption that the high-dimensional data lie on a low dimensional manifold $$ \mathcal{M}$$ [2–4]. In order to “learn” the manifold from a discrete set of data points, a graph is typically defined, where the graph nodes are the data points and the edges are determined according to some similarity notion. Since the manifold information is entirely captured by its Laplacian, the discrete counterpart, the graph Laplacian is used to build a low-dimensional embedding that respects the manifold in some sense [36]. To this end, common practice is to compute and exploit the spectral decomposition of the graph Laplacian. Diffusion maps is one of these methods, which constructs a random walk on the graph and represents the data points in a low-dimensional space preserving the neighborhood information [4].

Consider a set of data points $$ {\left\{y_i\right\}}_{i=1}^b$$, where $$ y_i\in {\mathbb{R}}^d$$ for i = 1, …, b. An undirected weighted graph $$ \mathcal{G}=(\mathcal{V},\mathcal{E},\mathit{W})$$ is constructed from the data points, where the vertex set is $$ \mathcal{V}=(y_1,y_2,...,y_b)$$ and the weights of the edges connecting two vertices are determined by a measure of similarity between any two data points, e.g., by

Next, a random walk P on the data points is constructed by normalizing the weight matrix W

Since P is similar to a symmetric and positive-define matrix, P has a biorthogonal right- and left-eigenvectors $$ {\{{\mathit{\phi }}_i,{\mathit{\upsilon }}_i\}}_{i=1}^b$$ with the eigenvalues 1 = λ₁ ≥ λ₂ ≥ … ≥ λ_b ≥ 0. Consequently, the spectral decomposition of P^t is given by

The diffusion distance $$ D_t^2(i,j)$$ between two data points y_i and y_j in the data set is defined by

The diffusion maps is defined by [4]

We remark that in many cases, due to the typical fast decay of the eigenvalues of P^t, l can be set to be smaller than d, thereby achieving dimension reduction. In addition, φ₁ is a constant vector and therefore is not used in the mapping.

It can be shown that the diffusion distance can be approximated by the eigenvalues and eigenvectors by [4]

Determining the dimension of data

A common problem in diffusion maps setting is how to choose the dimension l. The authors in [6] proposed Jackstraw to identify the number of principal components (PCs) in the context of principal component analysis (PCA) [37]. We present here a variant of Jackstraw, adapting it to diffusion maps.

Given a random walk P constructed from a set of b data points, the associated eigenvalues are 1 = λ₁ ≥ λ₂ ≥ … ≥ λ_b with corresponding right-eigenvectors $$ {\left\{{\mathit{\phi }}_i\right\}}_{i=1}^b$$. Collect the eigenvalues into a vector, denoted by λ, where λ = (λ₁, λ₂, …, λ_b)^⊤. Let $$ {\mathit{P}}_k^{*}$$ consist of the random permutation of P. Apply eigenvalue decomposition to $$ {\mathit{P}}_k^{*}$$ and obtain the corresponding eigenvalues vector $$ {\mathbf{\lambda }}_k^{*}$$. Repeat this shuffling procedure s times and obtain a set of vectors $$ {\left\{{\mathbf{\lambda }}_k^{*}\right\}}_{k=1}^s$$. The dimension of the representations is determined by

Note that the absolute values of the eigenvalues of $$ {\mathit{P}}_k^{*}$$ are considered because $$ {\mathit{P}}_k^{*}$$ is not necessarily symmetric and therefore its eigenvalues $$ {\mathbf{\lambda }}_k^{*}$$ are not guaranteed to be real.

Related work

Heat kernel signature

There are several shape analysis signatures obtained by spectral methods with different geometric properties such as isometry and deformation invariance [16, 17, 38, 39], related to the proposed method. One of the notable shape signatures is based on the heat diffusion on a shape, called Heat Kernel Signature (HKS) [16]. Broadly, the HKS is obtained by the eigenvalue decomposition of the heat kernel defined on the shape. In the context of our problem, since the heat kernel and the Laplace-Beltrami operator Δ_i share the same eigenbasis, and since the discrete graph Laplacian converges (point-wise) to the Laplace-Beltrami [40]

then, a discrete counterpart of HKS is given by

where

λ_k and ψ_k are the k-th eigenvalue and k-th eigenvector of the random walk, respectively, and t is the number of random walk steps on the graph. For more details on HKS, we refer the readers to [16].

Similarly to the DKS in Eq (12), the HKS can also be viewed as a low-pass filter. Observe that, for small t, the HKS approximates the DKS by Taylor expansion; for other t values, the weights assigned by the HKS decay faster than the weights of DKS, and therefore, the DKS gives more attention to finer structures.

Wave kernel signature

Another related shape signature is built by the wave function to the Schrödinger equation describing the quantum mechanical particles, called Wave Kernel Signature (WKS) [17]. Similarly to the HKS, the WKS is given by

where

with

λ_k and ψ_k are the k-th eigenvalue and k-th eigenvector of the random walk, and t is number of random walk steps on the graph. While the DKS and HKS are viewed as lowpass filters, we note that the WKS can be viewed as a band-pass filter. For more details, we refer the readers to [17].

Localization toy problem

To illustrate the challenge in the problem setting and the generality of the proposed solution, we present three simulations of different localization problems. The Matlab code is available in S1 Matlab Code.

Simulation 1

Consider 800 objects on a 2-sphere in $$ {\mathbb{R}}^3$$ that can be located at four different regions. Each region consists of 200 objects. The positions of the objects are measured by 5 sensors, giving rise to the following set of observations $$ {\left\{f_j\right(x_i\left)\right\}}_{j=1}^5\in {\mathbb{R}}^{100}$$, where j is the index of the sensor and i is the index of the object (position). Each sensor measures the position in d = 100 coordinates in the following way

where the standard deviation of the measurement depends on the distance between the position of the sensor and the position of the object $$ {\sigma }_i=20exp(-\parallel x_i-s_j{\parallel }_2^2)$$, I₁₀₀ is the identity matrix of size 100 × 100, ‖ ⋅ ‖₂ denotes the Euclidean norm and s_j denotes the 3D position of sensor j ∈ {1, …, 5}. In other words, each object position is captured by d = 100 realizations of a Gaussian random variable with variance that is proportional to the distance between the sensor position and the object position. Note that the positions are captured by the sensor through the variance, therefore, they are difficult to infer directly by the multi-feature observations.

In Fig 4A, we present our setting consisting of objects and sensors located on and near a sphere, respectively. The objects positions are marked by dots, the different regions are marked by different colors (red, black, blue and yellow), and the sensors positions are marked by green stars. At the bottom, we present the (high-dimensional) sensor observations. Each block consists of 100 × 200 scalar observations corresponding to the observations of a single sensor from each region, where d = 100 is the dimension of each observation and 200 is the number of the positions per region. Visually, it is evident that distinguishing between the different regions merely based on these observations is non trivial.

Fig 4

Illustration of localization toy problems.

(A) The sensors and objects locations on a sphere are marked by green stars and dots, respectively. The multi-sensor observations correspond to Simulation 1 (left), Simulation 2 (middle), and Simulation 3 (right). (B) Results of the application of our approach to Simulation 1: the SSDs obtained by the proposed method (left), the diffusion maps embedding (middle), and the localization confusion matrix obtained by a 10-fold cross-validation with an RBF SVM classifier (right). (C) A comparison between the localization accuracy obtained by the proposed method based on SSDs and the localization accuracy obtained based on the output of each sensor as well as the concatenation of the output from all the sensors. The localization accuracy is the number of correctly identified positions divided by the true number of total positions in each region.

For illustration purposes, we view the problem as a classification problem, where given the high-dimensional multi-feature observations, the task is to identify in which region the object at x_i resides.

The observations from each sensor consisting of 800 object positions from four regions were processed in two stages: first, a 3D embedding is constructed by applying diffusion maps to the high-dimensional observations, and second, an RBF SVM is applied to the obtained embedding in order to classify the region. To evaluate the classifiers, we perform a 10-fold cross-validation. A similar two-step procedure is applied to the concatenation of the observations from all the sensors.

Table 2 presents the resulting classification accuracy, i.e., the number of correctly classified positions divided by the total number of positions in each region. The presented results are obtained using a 10-fold cross-validation. We observe that none of the sensors enables an accurate classification. Moreover, we show in Table 2 that a naïve concatenation of the observations from all the sensors do not yield a good classification either.

Table 2

Localization accuracy from measurements of Simulation 1.

Region	Concatenated sensors	Sensor 1	Sensor 2	Sensor 3	Sensor 4	Sensor 5
Blue	43%	46%	18%	44.5%	58%	40.5%
Red	54%	15.5%	52.5%	45.5%	48.5%	55%
Black	46.5%	49.5%	29%	40%	37.5%	42%
Yellow	38.5%	33.5%	52%	30.5%	42.5%	50.5%

The performances obtained by the 10-fold cross-validations are presented. The percentages indicate the number of correctly classified positions divided by the true number of total positions in each class.

Seemingly, in order to mitigate the problem, we could simply represent the objects positions by a vector of the variance of the observations. However, it would require prior knowledge about the sensing model, whereas our approach is model-free.

In Fig 4B, we present the classification results obtained by the proposed method. In the diffusion maps embedding constructed from the SSDs, we observe a clear separation between the four regions. Finally, in the confusion matrix of the classification, we observe that the proposed method leads to significantly better classification results compared to the results in Table 2. This demonstrates the importance of taking into account the mutual-relationships between the sensor observations, rather than processing the nominal values of the observations directly, which in this case, give rise to correct identification of the four regions.

Simulation 2

The main purpose of this simulation is to demonstrate Proposition 2 from Binary hypothesis testing.

Consider n = 400 positions, $$ {\left\{x_i\right\}}_{i=1}^{400}$$, such that $$ x_i\in S^2\subset {\mathbb{R}}^3$$, which are sampled from two different regions on the sphere. Suppose that the positions of the objects follow a bimodal distribution as depicted in the middle-top figure in Fig 4A: an object is located in the blue region with probability $$ \frac{1}{2}$$ and in the red region with probability $$ \frac{1}{2}$$. The two regions represent the two hypotheses, $$ {\mathcal{H}}_1$$ and $$ {\mathcal{H}}_2$$, where each region consists of 200 positions. Note the symmetry in this setting, that is, the distances from the blue and red regions to the sensors are approximately the same.

Here, we have 12 sensor observations $$ {\left\{f_j\right(x_i\left)\right\}}_{j=1}^{12}$$, measuring the positions of the objects, located at s_j. The multi-feature observations are random samples from

where $$ f_j\left(x_i\right)\in {\mathbb{R}}^{100}$$, the standard deviation is $$ {\sigma }_j^i=20exp(-\parallel x_i-s_j{\parallel }_2^2)$$, and the covariance $$ {\mathbf{\Sigma }}_{k,l}\sim \left|\mathcal{N}(0,0.5)\right|$$ for $$ \lfloor \frac{k}{100}\rfloor ,\lfloor \frac{l}{100}\rfloor \in [1,\dots ,12]$$.

In this simulation, a direct computation can show that Assumptions (A.1) and (A.2) hold. Specifically, we note that

and

where Ω_r and Ω_b denote the sets of red or blue object locations, respectively. In addition, the conditions of Special Case 1 are satisfied, and thus, the total variation of these sensor observations is zero. In other words, using a single sensor is insufficient to distinguish between the two regions. Conversely, we show that since the SSDs take into account the covariance information between the sensors, they allow us to make this distinction.

In Fig 4C, we present a comparison between the classification results obtained by our approach using SSDs and the classification results obtained by using the output of each sensor as well as the concatenation of the output from all the sensors. The results are evaluated with a 10-fold cross-validation. We observe that the classification obtained by proposed algorithm is significantly better than the classification obtained using the “raw” sensor outputs.