MIAMI: Mutual Information-based Analysis of Multiplex Imaging data.

MOTIVATION: Studying the interaction or co-expression of the proteins or markers in the tumor microenvironment (TME) of cancer subjects can be crucial in the assessment of risks, such as death or recurrence. In the conventional approach, the cells need to be declared positive or negative for a marker based on its intensity. For multiple markers, manual thresholds are required for each marker, which can become cumbersome. The performance of the subsequent analysis relies heavily on this step and thus suffers from subjectivity and lacks robustness.
RESULTS: We present a new method where different marker intensities are viewed as dependent random variables, and the mutual information (MI) between them is considered to be a metric of co-expression. Estimation of the joint density, as required in the traditional form of MI, becomes increasingly challenging as the number of markers increases. We consider an alternative formulation of MI which is conceptually similar but has an efficient estimation technique for which we develop a new generalization. With the proposed method, we analyzed a lung cancer dataset finding the co-expression of the markers, HLA-DR and CK to be associated with survival. We also analyzed a triple negative breast cancer dataset finding the co-expression of the immuno-regulatory proteins, PD1, PD-L1, Lag3 and IDO, to be associated with disease recurrence. We demonstrated the robustness of our method through different simulation studies. AVAILABILITY: The associated R package can be found here, https://github.com/sealx017/MIAMI. SUPPLEMENTARY INFORMATION: The Supplementary Material is attached.

1 Introduction

In recent years, multiplex tissue imaging (Bataille ) technologies like, imaging mass cytometry (Ali ), multiplex immunohistochemistry (mIHC) (Tan ) and multiplexed ion beam imaging (MIBI) (Angelo ) have become increasingly popular for probing single-cell spatial biology. The technologies help in understanding the biological mechanisms underlying cellular and protein interactions in a wide array of scientific contexts. MIBI (Ionpath Inc.) (Angelo ) platform and the mIHC platforms such as Vectra 3.0 (Akoya Biosciences) (Huang ) and Vectra Polaris (Akoya Biosciences) (Pollan ) produce images of similar structure. In particular, each image is two-dimensional, collected at cell- and nucleus-level resolution and proteins in the sample have been labeled with antibodies that attach to cell membranes. We will refer to the antibodies as markers in the article. Typically, mIHC images have 6–8 markers, whereas MIBI images can have 40 or more markers.

Many of the above markers are surface or phenotypic markers (Shipkova and Wieland, 2012; Zola ), which are primarily used for cell type identification. Additionally, there are functional markers (Ijsselsteijn ) such as HLA-DR (Saraiva ), PD1, PD-L1 (Alsaab ) and CD45RO (Lee ) that dictate or regulate important cell functions. Both surface and functional markers are quantified as continuous valued marker intensities. However, in the traditional method, the interaction or co-expression effects of the markers are studied by binarizing them. For every marker, a threshold is chosen to indicate whether a cell (or, a pixel) in the tumor microenvironment (TME) (Binnewies ) is positive or negative for that marker. If a cell is positive for two markers, it implies that the markers have co-expressed or co-occurred in that cell. Next, for all the subjects considering every unique pair of markers, the proportion of cells positive for both the markers, the proportion of cells positive for only the first marker and the proportion of cells positive for only the second marker are computed. These proportions can then be used in hierarchical clustering (Murtagh and Legendre, 2014) to group the subjects. It can then be tested if the cluster labels correlate with clinical outcomes, such as disease recurrence or time to death (Jackson ; Johnson ; Koguchi ). The step of binarizing the marker expression profiles can be performed either using compatible commercially available software or manual assessment of the segmented images. For example, Johnson used AQUAnalysis® software (Dolled-Filhart ; McCabe ) for binarizing the functional markers, PD1, PD-L1, HLA-DR and IDO, and discovered that their co-expression predicted improved outcomes of anti-PD1 therapies in metastatic melanoma. On the other hand, Patwa used their own clinical expertise and careful evaluation of the segmented images to determine the binarizing thresholds and discovered the interaction between the markers, PD1, PD-L1, IDO and Lag3 to be associated with disease recurrence. It should be pointed out that such expression-thresholding or binarization is also popular in cell–cell communication or interaction analysis (Jin ) in the context of single-cell RNA sequencing data. As an example, Armingol defined a pair of cells to be interacting if the expressions of both ligand and receptor (Wong ) in those cells exceed certain chosen thresholds.

The manual threshold selection process for multiple markers can be extremely challenging and is subjective. Alternatively, commercially available softwares, such as AQUAnalysis® and inForm (Dolled-Filhart ; Kramer ), can be used for automated threshold selection in the context of a few particular data types. However, these softwares generally follow a ‘black box approach’ and can be difficult to interpret. On top of these difficulties, many authors have criticized binarizing continuous random variables in general due to the resulting loss of power (Altman and Royston, 2006; Irwin and McClelland, 2003; Seal ).

In this article, we propose a threshold-free approach for studying marker co-expression. We treat the marker intensities as continuous random variables and use their marginal and joint probability density functions (PDF’s) to construct a metric of co-expression based on mutual information (MI) (Cover and Thomas, 2006). Unlike the correlation coefficient, MI is capable of capturing non-linear patterns of dependence between the markers and is easily extendable for more than two markers. However, as the number of markers increases, computing the joint PDF becomes increasingly challenging which makes the computation of MI infeasible as well. Therefore, we use a slightly different formulation of MI known as Euclidean quadratic mutual information (EQMI) (Principe ) which has a similar interpretation but can be computed more efficiently. The computation algorithm is discussed in Principe (2010) with a simpler assumption that we further generalize. The vector of estimated values of the EQMI of all the subjects is tested for association with clinical outcomes. With the proposed method, we analyzed an mIHC lung cancer dataset (Seal ) finding that a higher co-expression of the markers, HLA-DR and CK was significantly associated with better 5-year overall survival. We analyzed a MIBI triple-negative breast cancer (TNBC) dataset (Keren ) studying the co-expression of two sets of functional markers, (a) HLA-DR, CD45RO, H3K27me3, H3K9ac and HLA-Class-1, and (b) PD1, PD-L1, Lag3 and IDO, which are also known as immuno-regulatory proteins (IRP’s). We found the co-expression of the IRP’s to be significantly associated with recurrence. We demonstrated the robustness of our method over the existing approaches through different simulation studies.

2 Materials and methods

Suppose there are p markers and N subjects with the j-th subject having n cells. Let X denote the expression of the k-th marker in the i-th cell of j-th subject for , and Note that we focus on cell-level data in this article but the framework is readily usable on pixel-level data as well. Let denote a subject-level outcome vector and C be an N × S matrix of subject-level covariates. Next, we discuss the existing and proposed methods and also summarize them in Figure 1.

Fig. 1.

Comparison of the workflow of the proposed method with the traditional method. We used segmented cell-level data in this article but the method is applicable on a pixel-level data as well

2.1 Traditional thresholding-based approach to study co-expression of the markers

For every marker k, we choose a cutoff t to define cell i of subject j as positive for that marker if The choice of t can be guided by prior biological insight or by careful inspection of the marker intensity profile. For example, extreme quantiles like (e.g. the 90th or 95th percentiles) can serve as viable thresholds, as we see in the simulations. However, in most real datasets, it requires user-defined thresholds for achieving a meaningful and interpretable conclusion. A cell can be positive for multiple markers. If a cell is positive for a pair of markers (k, k), it would imply that these markers have co-expressed in that particular cell. For every subject, compute the proportion of such double-positive cells, denoted by , the proportion of cells positive for only the first marker, denoted by , and the proportion of cells positive for only the second marker, denoted by , for and Next, the subjects are classified into two or more groups based on their vectors of proportions using hierarchical clustering. Note that if it is biologically relevant, instead of a pairwise analysis, one can also study multiple markers jointly, e.g. with four markers, one can count the cells which are either or any of the possible combinations and group the subjects based on these proportions.

Suppose that the subjects are grouped into M clusters. Let be an N × M matrix of the cluster labels. Z is a vector corresponding to the subject j with Z = 1 if the j-th subject belongs to group m and 0 otherwise. In most common practices, M = 2 is considered. When Y is a continuous outcome, a standard multiple linear regression model with Z as a predictor can be written as where are fixed effects and ϵ is an error vector following multivariate normal distribution (MVN) with mean 0 and identity covariance matrix The null hypothesis, , can be tested using the Wald test or likelihood ratio test (LRT) (Gourieroux ). Similarly, when Y is a categorical outcome, a logistic or multinomial logistic regression model (Kwak and Clayton-Matthews, 2002) can be considered.

Next, we consider the case of Y being a right-censored failure time outcome. Let the outcome of the j-th individual be , where T is the time to event and U is the censoring time. Let be the corresponding censoring indicator. Assuming that T and U are conditionally independent given the covariates for , the hazard function for the Cox proportional hazards (PH) model (Andersen and Gill, 1982) with fixed effects can be written as, where is the hazard of the j-th subject at time t, given the vector of covariates C and the cluster label Z and is an unspecified baseline hazard at time t. To test the null hypothesis: , an LRT (Therneau, 1997) can be considered.

2.2 Proposed method: mutual information-based analysis of marker co-expression

2.2.1 Theory of mutual information

MI is an information theoretic measure of dependence between two or more random variables (r.v.’s). In contrast to the linear correlation coefficient, it captures dependences which do not manifest themselves in the covariance (Kraskov ). For two random variables, R1 and R2 with sample space , marginal PDFs, f1, f2 and joint PDF f12, the MI is defined as the following:

MI has been used as a tool for feature selection in many different contexts (Hoque ; Liu ; Song ; Yang and Moody, 1999). The measure can be easily generalized for more than two random variables. However, due to the curse of dimensionality (Langrené and Warin, 2019), it becomes extremely difficult to estimate the joint PDF and hence, the MI, as the number of random variables increases.

Xu (1998) and Principe looked into alternative definitions of MI that would remain conceptually similar but can be computed efficiently. They discussed two measures, EQMI and Cauchy–Schwarz quadratic mutual information (CSQMI), defined as below,

It is trivial to verify that with equality occurring if and only if R1, R2 are independent, i.e. for any r1, r2. Using the Cauchy–Schwarz inequality, we have with equality happening if and only if R1, R2 are independent. Xu (1998) argued that EQMI shares more properties, such as convexity with respect to the PDF’s, of the traditional MI compared to CSQMI and thus, is better suited as a dependence measure. It should also be noticed that the form of EQMI is very similar to the distance covariance measure proposed by Székely . There is one more generalized measure of dependence, known as kernel canonical correlation analysis (Huang ), which does not share forumlaic similarity with EQMI but can potentially serve as an alternative for detecting non-linear dependence patterns. In this article, we focus on EQMI and propose its usage in the co-expression analysis of the markers in general multiplex imaging datasets used in the study of spatial biology of TME.

2.2.2 Formulation of EQMI

For every subject j, we assume that the expression of marker k is a continuous random variable, denoted by X, with sample space . X is observed in n cells as, . Suppose there are p = 2 markers and for every subject j, denote their joint PDF as and their marginal PDFs as and , respectively. Following Equation (2), the EQMI between the markers 1 and 2 for subject j can be defined as, where , , and . can be interpreted as a generalized measure of co-expression of the markers, capable of capturing non-linear dependences. A large value of will imply that the markers 1, 2 have significantly co-expressed in the TME of subject j. is bounded below by 0 for every j but there is no common upper bound for different j’s. Therefore, to compare ’s across different subjects, we need to appropriately standardize the values so that they lie in the same scale. We define a new measure as, We observe that lies between because , and has a similar interpretation as .

It is intuitive to generalize the measure for any markers. For subject j, letting to be the joint PDF and to be the marginal PDFs, can be defined as,

Notice that to estimate when the true PDFs are unknown, a naive approach would be to estimate the PDFs first by using a kernel density estimation approach (Silverman, 1981). Then, we use the estimated PDFs to compute the terms V, V and V via numerical integration (Davis and Rabinowitz, 2007). However, such an approach would be computationally infeasible for a large p and will defeat the purpose of considering EQMI instead of the standard form of MI.

can be estimated efficiently as, , where and . G stands for a Gaussian kernel with bandwidth parameter h. This clever way of estimating the terms, and is described in Principe (2010) with a simpler assumption that hs are equal for all k, i.e. h = h for We provide the general derivation (i.e. , for ) and other associated details in the Supplementary Material. For choosing the optimal values of hs, we generally consider the diagonal multivariate plug-in bandwidth selection procedure (Chacón and Duong, 2010; Wand ). However, when p is large (p > 6), to avoid computational deadlock we suggest using Silverman’s rule of thumb (Silverman, 1981) for choosing hs individually.

2.2.3 Using mutual information in association analysis

Denote the vector of estimated values of as , where . Our goal is to test if E is associated with the clinical outcome Y. When Y is a continuous outcome, a standard multiple linear regression model with E as a predictor can be written as, where are fixed effects and ϵ is an error vector following MVN with mean 0 and identity covariance matrix After estimating the parameters, the null hypothesis, , can be tested using the Wald test. Note that we can also add higher order terms of E, such as as in a polynomial regression model (Ostertagová, 2012) to account for non-linear relationship between Y and E. Similarly, when Y is a categorical outcome, a logistic or multinomial logistic regression model (Kwak and Clayton-Matthews, 2002) can be considered.

Next, we consider the case of Y being a survival or recurrence outcome. Using the same definitions and conditional independence assumptions of T, U and covariates as in Section 2.1, the hazard function for the Cox PH model can be written as, where is the hazard of the j-th subject at time t, given the vector of covariates C and the predictor E and is an unspecified baseline hazard at time t. To test the null hypothesis: , an LRT can be considered.

In a two-marker scenario (i.e. p = 2), one can also treat the absolute value of the Pearson correlation as a measure of marker co-expression and use it as E in the earlier equations for testing association with the outcome. However, as we later demonstrate, using correlation instead of EQMI can be sub-optimal in many cases. It is also difficult to generalize to more than two markers.

3 Real-data analysis

We applied our method to two real datasets, an mIHC lung cancer dataset (Seal ) from Vectra 3.0 platform and a MIBI TNBC dataset (Keren ). We also applied the traditional thresholding-based method on both the datasets. Since it was hard to decide the optimal thresholds for binarizing the markers, we ran the method for varying values of the thresholds. For every marker, concatenating the intensity data of all the subjects, we computed the median, 95% and 99% quantiles. Next, three different thresholding-based methods using these quantiles were considered, respectively referred to as Median-Thresholding, Threshold 1 and Threshold 2. Note that Threshold 1 and 2 both captured the difference in the tails of the distributions, whereas Median-Thresholding captured the difference in the centers. For the first dataset, we also performed the correlation-based association analysis, referred to as Corr.

3.1 Application to mIHC lung cancer data

In the lung cancer dataset, there were 153 subjects each with 3–5 images (in total, 761 images). Two subject-level covariates, age and sex were available. For every subject, the images were non-overlapping and from the same TME region. After segmenting the images, the subjects had varying number of identified cells (from 3755 to 16 949). We worked directly with the cell-level data as described next. The cells come from two different tissue regions: tumor and stroma. They were classified into either of the six different cell types: CD14+, CD19+, CD4+, CD8+, CK+ and Other, based on the expression of the phenotypic markers, CD19, CD3, CK, CD8 and CD14. These markers (and the cell-types) usually have clinical meaning, e.g. CK is a type of a tumor cell marker and CD4 and CD8 are T-cell or cytotoxic T-cell markers. Apart from these, a functional marker HLA-DR (also known as MHCII), was measured in each of the cells. Johnson classified the subjects into two groups, (a) MHCII: High and (b) MHCII: Low based on the proportion of tumor cells that are also positive for HLA-DR (i.e. cells). They discovered that group (a) had significantly higher 5-year overall rate of survival (reported P-value of 0.046).

Note that having a large number of CK+HLADR+ tumor cells implies that these two markers had co-expressed in a lot of the tumor cells. In light of that, we studied if the degree of co-expression of these markers in the tumor cells, as quantified by our method, was associated with the survival. Considering the tumor cells of every subject j, we first estimated between the two markers, HLA-DR and CK. Next, we tested the association of E with 5-year overall survival using the Cox-PH model from Equation (5). The coefficient γ was −7.26 with the P-value of the LRT being 0.0286. Thus, subjects with high co-expression of the markers in the tumor cells were more likely to survive. The result was thus consistent with Johnson ’s finding. The estimated coefficient, hazard ratio (HR) and P-value of all the methods are listed in Table 1 and further details such as confidence interval of the HR, are provided in the Supplementary Material. The correlation-based analysis (Corr) yielded a negative coefficient estimate but was not statistically significant (at level 0.05). Out of the thresholding-based methods, only Threshold 2 had a significant P-value. It demonstrated that the traditional thresholding-based method could vary heavily based on the choice of the thresholds leading to utterly different conclusions. On the other hand, inference using our method would be robust.

Table 1.

Estimated coefficient, hazard ratio (HR) and LRT P-value for testing association with 5-year overall survival using different methods in the mIHC lung cancer dataset

Method	Coefficient (HR)	LRT P-value
EQMI^*	−7.26 (0.0007)	0.0286
Corr	−1.38 (0.2512)	0.1838
Median-thresholding	−0.48 (0.6189)	0.0495
Threshold 1	−0.58 (0.5585)	0.0575
Threshold 2	−1.03 (0.3555)	0.0098

We should point out that CK was a phenotypic marker, and studying its co-expression with a functional marker HLA-DR might not have much clinical relevance. However, the goal was to demonstrate that our method obviates the need to binarize the continuous-valued marker expression profiles and is potentially applicable even when one of the markers is a phenotypic marker.

3.2 Application to TNBC MIBI data

The TNBC MIBI dataset (Keren ) had 38 subjects, each with one image. There were 201 656 cells in total with each of the images having varying numbers of cells (between 1217 and 8212). There were 49 markers in total, the majority of which were lineage or phenotypic markers, used primarily for cell-type identification. We were interested in two sets of functional markers, (a) HLA-DR, CD45RO, H3K27me3, H3K9ac and HLA-Class-1, and (b) PD1, PD-L1, Lag3 and IDO, also known as IRP’s. The reason for concentrating on these two sets of markers was the findings of Patwa . Employing the thresholding-based method with a set of very carefully chosen thresholds, they concluded that the pair-wise co-expression of the functional markers from the sets (a) and (b) were negatively associated with two clinical outcomes, namely disease recurrence and time to death (survival). For set (a), they did not find any statistical significance for either of the clinical outcomes (at level 0.05), whereas, for set (b), they were able to find statistical significance in the association test with recurrence (reported P-value of 0.0058).

For the co-expression analysis with the five markers from the set (a), we looked into all possible () two-way and higher-order combinations of the markers. There were four different types of combination, namely pair (10), triplet (10), quadruplet (5) and quintuple (1). We computed for each combination of the markers and tested for association with recurrence and survival. In Table 2, we list the results for five marker combinations for which the lowest P-values were observed. We noticed that at a level of 0.05, several marker-combinations were found to be associated with recurrence. All the estimated coefficients were negative, implying that higher co-expression of the markers decreased the chance of disease recurrence. However, once we adjusted the P-values for multiple testing correction using Bonferroni’s method (Bonferroni, 1936), only the combinations ‘HLA-DR, CD45RO’ and ‘HLA-DR, CD45RO, H3K9ac’ remained significant. Note that, we compared the P-values of every type of combination separately, meaning that for marker pairs and triplets, we compared the P-values at level 0.05/10 since the numbers of pairs and triplets were both 10, and for quadruplets, at level 0.05/5 since the number of quadruplets was 5. The marker-combinations were not independent and thus, a Bonferroni correction probably has been overly conservative in this case. For the survival outcome, we did not detect any statistical significance. But, the negative coefficient estimates hinted at a possible association of better rate of survival with higher co-expression.

Table 2.

Estimated coefficient and LRT P-value for testing association with recurrence and survival for five combinations of the markers from sets (a) and (b) with the lowest P-values, obtained by the proposed method in the MIBI dataset

Clinical	Marker	Coefficient	LRT
Outcome	Combination		P-value
Recurrence	HLA-DR, CD45RO	−32.97	0.0022
	HLA-DR, CD45RO, H3K9ac	−12.48	0.0051
	HLA-DR, CD45RO, H3K27me3	−6.73	0.0245
	CD45RO, H3K9ac	−41.53	0.0329
	HLA-DR, CD45RO, HLA-Class-1	−7.20	0.0421
Survival	HLA-Class-1, H3K27me3	−12.75	0.1010
	HLA-DR, H3K9ac	−8.27	0.1542
	HLA-DR, CD45RO	−8.35	0.1583
	HLA-DR, HLA-Class-1	−29.69	0.1630
	HLA-DR, CD45RO, H3K9ac	−0.74	0.3158
Recurrence	PD1, PD-L1	−9.61e + 03	0.0046
	PD1, PD-L1, IDO	−5.54e + 02	0.0065
	PD1, PD-L1, Lag3, IDO	−4.37e + 02	0.0069
	PD-L1, Lag3, IDO	−8.94e + 02	0.0084
	PD1, PD-L1, Lag3	−1.91e + 02	0.0090
Survival	PD-L1, Lag3	−1.20e + 02	0.0103
	Lag3, IDO	−1.14e + 02	0.0302
	PD-L1, Lag3, IDO	−5.67e + 02	0.0449
	PD-L1, IDO	−6.85e + 02	0.0490
	PD1, PD-L1, Lag3, IDO	−2.15e + 02	0.0586

For the co-expression analysis with the four IRP’s from the set (b), we looked into all possible () two-way and higher order combinations. There were four different types of combination, namely pair (6), triplet (4) and quadruplet (1). In Table 2, we list the results for five marker combinations for which the lowest P-values were observed. At a level of 0.05, many of the marker-combinations were found to be associated with both recurrence and survival. Again, all the estimated coefficients were negative, implying that higher co-expression of the markers decreased the chance of disease recurrence. Upon correcting the P-values for multiple testing using Bonferroni’s method, all of the marker combinations listed in Table 2 remained significant for recurrence while for survival, none of them remained significant. We compared the P-values of every type of combination separately, meaning that for marker pairs we compared the P-values at level 0.05/6 since the number of pairs was 6, for triplets we compared the P-values at level 0.05/4 since the number of triplets was 4, and for quadruplet, at level 0.05 since there was just a single quadruplet. Thus, we also arrived at a similar conclusion as Patwa that the inter-play or co-expression of the IRP’s were significantly associated with recurrence and possibly also with survival. One added novelty of our method was that one could easily pinpoint which of the combinations of the IRP’s had the most impact.

It should be kept in mind that the sample-size for this dataset was quite small with only 16 events for recurrence and 15 for survival. It might have affected the overall inference which was based on asymptotic distributional properties of the test statistics. For the same reason, we mainly focused on the sign of the coefficient estimates but not their CI’s in this particular case. We also applied the simple thresholding-based methods, Median-Thresholding, Threshold 1 and Threshold 2 as described earlier using both the sets of markers. We found only a single statistically significant result which was for Median-Thresholding using set (b) in the association test with recurrence. Associated tables are provided in the Supplementary Material.

4 Simulation

Next, we compared the performance of the -based association analysis with the correlation-based association analysis and the thresholding-based methods in different simulation setups. We assumed that there were two groups of subjects, one with subjects having high marker co-expression and the other with subjects having low or almost zero marker co-expression. In Section 4.1, we considered two markers sharing linear and non-linear patterns of co-expression relationship. In Section 4.2, we considered three markers sharing varying degree of linear co-expression relationship.

4.1 Simulation with two markers

4.1.1 Simulation using Gaussian copula

We replicated the characteristics of the lung cancer dataset in this simulation. The mean marginal distribution of the markers HLA-DR and CK across the subjects could be approximated by Beta distribution respectively with parameters, and (refer to the Supplementary Material). We used a Gaussian copula (Masarotto and Varin, 2012) to simulate correlated intensity data for two markers which had the above marginal Beta distributions. The simulation strategy was as follows,

A random number I between 0 and 1 was chosen with probability 0.5 each, respectively standing for group (1), whose subjects had high co-expression of the markers and group (2), whose subjects had mild to none co-expression of the markers. It assigned j-th subject to either of the two groups.

The intensity vector of two markers, for every individual j was simulated as follows,

If I = 0, simulate a correlation parameter ρ from , or else simulate ρ from .

Consider a correlation matrix, and simulate for

Compute , where denotes the cumulative distribution function of the standard normal distribution.

Perform inverse transformation as,

, where F1 and F2 are Beta distributions with parameters and . Refer to the Supplementary Material for plots of the true joint densities of the markers for two groups of subjects.

The clinical outcome of j-th subject is simulated as, , where This step is repeated 100 times to generate 100 different datasets having different Y vectors but the same intensity vectors, X1 and X2. All the methods are applied on these 100 datasets and empirical power is computed.

Steps were repeated 20 times and the mean empirical power of different methods were displayed for varying values of the number of cells and the number of individuals (N) in Figure 2. The -based association analysis (EQMI) and the correlation-based association analysis (Corr) achieved comparable performance in all the cases. This particular simulation strategy inherently assumed that the dependence between the markers was linear. Thus, the estimated values of the and the correlation shared an almost one-to-one relationship making the association analysis using either of them equivalent. Median-Tresholding and Threshold 1 showed similar performance, whereas Threshold 2 had consistently lower power. It showed the importance of choosing a proper threshold in the traditional thresholding-based method for achieving a reasonable performance. All the methods expectedly had the least power when N was the smallest, whereas the value of did not have any major impact. It implied if the co-expression pattern was well captured even through a smaller number of cells, most of the methods would perform well.

Fig. 2.

The figure displays the power of different methods under different simulation scenarios from Section 4.1 with two markers for varying numbers of subjects (N) and cells (). On the x-axis, the fixed effect size β was varied from low to high

4.1.2 Simulation with squared marker co-expression relationship

The last simulation strategy essentially assumed a linear pattern of co-expression between the markers. Next, we simulated a non-linear pattern of co-expression between the markers where we would expect the correlation-based association analysis (Corr) to perform worse since it could only capture a linear dependence pattern. Steps 1 and 3 of the last simulation strategy were kept the same here and the marker-intensity simulation step, i.e. step 2 was changed as follows,

2a. If I = 0, simulate from Unif, and e from Unif. Construct as,

2b. If I = 1, independently simulate both from Unif.

From Figure 2, we noticed that the -based association analysis performed the best in all the cases. Threshold 1 and 2 achieved comparable performance for large value of , whereas Median-Thresholding mostly yielded poor performance. Note that for both the groups of subjects, the correlation between the markers were close to zero as the dependence pattern was non-linear, squared to be specific. Expectedly, the correlation-based association analysis (Corr) had almost no power in every case. The simulation strategy demonstrated why using a generalized measure of co-expression such as would be more optimal in many cases.

4.1.3 Simulation with circular marker co-expression relationship

In the last simulation setup, the thresholding-based methods performed well despite the marker co-expression pattern being non-linear. Next, we looked into a slightly more complicated co-expression relationship for which the thresholding-based approach would suffer. Steps 1 and 3 of the last two simulation strategies were kept the same here and the Step 2 was changed as follows,

2a. If I = 0, simulate from Unif, e from Unif and a random number s between −1 and 1. Construct as,

2b. If I = 1, independently simulate both from Unif.

From Figure 2, we noticed that the -based association analysis performed the best in all the cases. The correlation-based association analysis (Corr) expectedly performed the worst. Threshold 1, unlike the last simulation, performed significantly worse. In this simulation setup, the subjects of one group had a circular pattern of marker co-expression and the others had almost zero co-expression. Recall that in a two-marker scenario, the thresholding-based methods depended on computing the proportions of the cells positive for both the markers and of the cells positive for only one of the markers (Section 2.1). The difference between these proportions across the two groups of subjects became negligible under this setup which made it difficult distinguishing between them. This explained the overall poor performance of all the thresholding-based methods.

4.2 Simulation with three markers

Next, we considered three markers and simulated varied degree of linear dependence between them using Gaussian copula. We only performed the -based association analysis and the thresholding-based methods in this case. There were again two groups of subjects respectively with high and low co-expression. The simulation strategy was as follows,

A random number I between 0 and 1 was chosen with probability 0.5 each, respectively standing for groups (1) and (2).

The intensity vector of three markers, for every individual j was simulated as follows,

Consider a correlation matrix, . For the subjects in group (1), the off-diagonal elements of the correlation matrix would have high values, whereas they would have low values for the subjects in group (2). We considered three different cases with varying differences between the correlation matrices of the two groups.

Case (a): If I = 0, the correlation parameters and were independently simulated from . Otherwise, they were independently simulated from .

Case (b): Regardless of the value of I, was kept to be 0. If I = 0, and were independently simulated from . Otherwise, these two were independently simulated from .

Case (c): Regardless of the value of I, and were kept to be 0. If I = 0, was simulated from . Otherwise, it was simulated from .

Simulate and compute .

Perform inverse transformation as, , where F1, F2 and F3 are the Beta distributions respectively with parameters and .

The clinical outcome of j-th subject was simulated as, , where This step was repeated 100 times to generate 100 different datasets having different Y vectors but the same intensity data, X1, X2 and X3. All the methods were applied on these 100 datasets and empirical power was computed.

Steps 1–3 were repeated 20 times and in Figure 3, the mean empirical power of the methods were displayed. The power of the methods were quite low when N was small. The -based method outperformed both the thresholding-based methods, Threshold 1 and 2 in every case. Threshold 2 had little to no power in most of the cases. Note that, the cases (a), (b) and (c) differed in how different the marker co-expression pattern of the two groups were. The difference between the marker co-expression pattern of the two groups was the largest in case (a) since all the three correlation parameters, and were different across the groups. The difference was the smallest in case (c) as two of the three correlation parameters, and were kept to be 0 in both the groups. Quite expectedly, the power of the methods decreased going from case (a) to case (c), as the difference between the groups of subjects reduced. The decrease was more prominent with Threshold 1 suggesting the method’s lack of robustness.

Fig. 3.

The figure displays the power of different methods under different cases from Section 4.2 with three markers for varying numbers of subjects (N) and cells (). On the x-axis, the fixed effect size β was varied from low to high

5 Discussion

In multiplex imaging data, studying the interaction or co-expression of multiple functional markers in the cells of the TME can be crucial for subject-specific assessment of risks. The traditional approach requires a complex step of binarizing the continuous valued marker expression profiles which is prone to subjectivity and can be sub-optimal in many scenarios. The complexity gets exacerbated as the number of markers increases. In this article, we propose a method for studying interaction or co-expression of multiple markers based on the theory of mutual information (MI). We treat the subject-specific intensity or expression of every marker as a continuous random variable. We determine how much the markers have co-expressed in the TME of a particular subject by computing a measure known as EQMI, comparing the estimated marginal and joint PDFs of the markers. The formula of EQMI has a similar interpretation as the standard formula of MI but allows a more efficient computation. We adopt and generalize an existing algorithm for computing EQMI that does not require explicitly estimating the joint PDF of the markers, a step which becomes increasingly intractable as the number of markers increases. Next, the subject-level EQMI values are tested for association with the clinical outcomes. The proposed method is free from the subjectivity bias of the traditional thresholding-based method and is readily applicable with any number of markers.

We applied the proposed method to two real datasets, one mIHC lung cancer dataset and one MIBI triple negative breast cancer dataset. In the former, we found high co-expression of the markers, HLA-DR and CK to be associated with the 5-year overall survival of the subjects. In the latter, we found high co-expression of the proteins, PD1, PD-L1, IDO and Lag3 to be associated with disease recurrence. We evaluated the performance of our method through several simulation scenarios with two and three markers. In the scenarios with two markers, we showed that all the methods perform well and close to each other if the pattern of dependence (co-expression) between the markers is linear. However, with a more complex non-linear dependence pattern, only the proposed method could achieve respectable power. In the scenarios with three markers, we found that the proposed method performed consistently better than the thresholding-based method and showed superior robustness.

As we have shown in the simulation studies, EQMI can capture both linear and non-linear patterns of co-expression between the markers very well. However, the measure is not well suited for capturing the differences between the patterns. For example, it may happen that one subject has a linear pattern of co-expression, whereas some other subject has a non-linear pattern. The EQMI for both the subjects can be very similar, making it hard to distinguish between them. As a part of our future direction, we would like to improve the method by detecting and incorporating the type of the co-expression pattern. With more than two markers, we studied the co-expression patterns of all possible combinations of the markers and declared significance based on P-values corrected by Bonferroni’s method. However, in future, we would like to explore the causal direction between the markers which can then be used to determine a smaller and optimal set of markers and would obviate the need of exploring all possible marker-combinations. In this article, we have not used any information on the spatial locations of the TME cells. As a future direction, we would like to study the MI between the spatial information and the marker expression profiles with a goal to detect spatially variable markers and their spatial patterns.

Our method is available as an R package named MIAMI at this link, https://github.com/sealx017/MIAMI. The package is readily applicable to any multiplex imaging dataset which has cell-level intensity data on two or more markers. In future, we would like to further augment the package’s capability by incorporating a pixel-level analysis as well.

Click here for additional data file.