Co-localization features for classification of tumors using mass spectrometry imaging

Statistical modeling of mass spectrometry imaging (MSI) data is a crucial component for the understanding of the molecular characteristics of cancerous tissues. Quantification of the abundances of metabolites or batch effect between multiple spectral acquisitions represents only a few of the challenges associated with this type of data analysis. Here we introduce a method based on ion co-localization features that allows the classification of whole tissue specimens using MSI data, which overcomes the possible batch effect issues and generates data-driven hypotheses on the underlying mechanisms associated with the different classes of analyzed samples.


INTRODUCTION
Mass spectrometry imaging (MSI) has been used to investigate the spatial distributions of molecular species within tissue specimens. In particular, it has shown promising results to serve as a technological platform for the study of chemical heterogeneity associated with the samples of interest. In the particular case of cancer tissue, established protocols exist that are tailored for the imaging analysis of human biopsy specimens (either freshly prepared -frozen tissue -or -formalin fixed paraffin embedded (FFPE), etc. - (Rubakhin et al., 2005;Patel, 2007)) using either desorption electrospray ionization (DESI) (Takáts, Wiseman and Cooks, 2005; Justin M. , matrix-assisted laser desorption mass spectrometry (MALDI-MS) (Hillenkamp et al., 1991;Spengler, 1997;Spengler and Hubert, 2002;Cornett et al., 2007;Lemaire et al., 2007;Casadonte and Caprioli, 2011) or secondary ionization mass spectrometry (SIMS) (Colliver et al., 1997;Chandra, 2004;Johansson, 2006;Debois et al., 2009) as the main predominant MSI analytical techniques. MSI spectra obtained in this way have been used to associate ionic species with histological features, to find biochemical differences between healthy or diseased (cancerous, necrotic, etc.) tissue regions or to study the tissue distribution of drugs (J. M. Nilsson et al., 2015;Karlsson and Hanrieder, 2017) and to relate small molecule metabolites and lipids with the expression of protein or other molecular markers (Rawlins et al., 2015). Although some of these methods involve elaborate sample preparationsuch as MALDI -, they also suffer from the usual analytical challenges of mass spectrometry-based analysis (ion suppression/matrix effects, intensity drifts, batch effects) (Taylor, Dexter and Bunch, 2018). Furthermore, the inherent complexity of solid tissue analysis brings extra difficulties regarding devising quality control materials and procedures for assessing ionization efficiency and matrix effects, platform stability and intensity drift. For that reason, the capacity to measure and correct the impact of batch and run-order effects is still underdeveloped compared to better established metabolomics and proteomic techniques (Gregori et al., 2012;Fages et al., 2013;Ellis, Bruinen and Heeren, 2014;Goh, Wang and Wong, 2017;Sánchez-Illana et al., 2018).
Although workflows which use peak intensities for classification or regression analyses are widespread in the literature, their effectiveness may be influenced by the magnitude of the batcheffects and other uncontrolled systematic variations. Without any means of the batch-effect correction, the relationship between the true molecular properties of the analyzed samples and the observed and molecular patterns becomes unreliable. Additionally, due to the nature of the MSI data, combining the information present in a large number of pixel spectra to discriminate entire tissue specimens can be challenging. For instance, studies have shown that MSI data can be used to distinguish different cancer types by selecting a random subset of pixel spectra from a region of interest (Meding et al., 2012), or to determine the risk of metastatic relapse if a certain number of pixel spectra are predicted according to determined labels (Abdelmoula et al., 2016).
A possible alternative approach, employed in the presented work, exploits more extensively the spatial information about the molecular distributions contained in MSI data, representing the molecular information of individual samples by their co-localization patterns. In this way, the whole sample can be expressed in its entirety by a single vector of co-localization features. Additionally, possible issues caused by batch-to-batch variability (Goh, Wang, and Wong, 2017;Bemis et al., 2018) can be mitigated, since co-localization is not affected by systematic variations of ion intensities across multiple samples that can be caused by instrumental or ion source/ion transfer conditions. In general, the biological rationale is that ion co-localization may reflect the underlying biochemical processes that characterize the various classes of samples of interest. A very similar idea has been employed for molecular imaging studies at the tissue and cellular scale, for instance in the case of fluorescence microscopy (Bermudez-Hernandez et al., 2017;Lamberts et al., 2017). Conceptually, the MSI data is converted into a graph, and then network features are used for supervised classification, in an analogous fashion to (Amoroso et al., 2018), where the complex networks framework is employed for classifying patients with neurodegenerative disease. Specifically, Pearson's correlation has been used as a suitable measure for determining the degree of co-localization, in particular in fluorescence microscopy (Manders et al., 1992).
In this paper, we show how the co-localization features can be used to classify three types of cancer specimens, and also provide insight into the possible mechanisms underlying their biochemical differences.

MATERIALS AND METHODS
Sixty-nine specimens from three distinct cohorts of human tumors were cryosectioned and subsequently analyzed by DESI-MSI. The full dataset consisted of twenty-eight breast tumors, part of the dataset published in (Guenther et al., 2015), sixteen colorectal adenocarcinomas, part of the dataset published in (Inglese, Strittmatter, et al., 2018) and twenty-five ovarian cancer samples, part the dataset published in (Dória et al., 2016). All the raw spectra were acquired in the negative ion mode using a custom DESI sprayer ion source interfaced to a high-resolution orbital trapping mass spectrometer [Exactive, Thermo Fisher Scientific] in the range of 150 -1,000 m/z using the described experimental and instrumental conditions. The mass spectrometry images were therefore split into two random subsets, denoted cross-validation set and test set consisting of sixteen samples per tumor type, and twelve breast plus nine ovarian cancer samples, respectively.

Pre-processing
All the raw spectral data were first converted into imzML format (centroided mode) (Schramm et al., 2012) and their m/z values were re-calibrated using a LOESS model fitted on the The re-calibrated image samples were independently pre-processed through a workflow consisting of a) TIC scaling intensity normalization, b) peak matching, c) log-transformation using MALDIquant for R (Gibb and Strimmer, 2012;Gibb, Franceschi and Gibb, 2017). The choice of TIC scaling normalization was based on the fact that correlations are insensitive to the domain where the measure is defined. This normalization method allowed us to consider the colocalization of the relative abundances of the ions, reducing the possible effects of systematic variations associated with more quantitative approaches. Peak matching was performed in a way so that no sample could have contributed to the matched m/z with more than one peak. No intersample normalization was applied to this dataset.
Tissue-unrelated peaks were filtered using the global reference and pixel count filters of the SPUTNIK package for R (Inglese, Correia, et al., 2018). The parameters used are reported in Table 1.
To match the peaks from the different images, peak matching was applied to the average of m/z values corresponding to the non-zero intensity peaks across all the pixels spectra of the cross-validation set. Only the peaks detected in all the images were retained for the statistical modeling (Table 2). This strict constraint was necessary to ensure that the different co-localization patterns among the tumor types were involving the molecules observed in all the samples, strengthening their biological significance. The final common m/z vector was de-isotoped using the procedure described in (Inglese, Strittmatter, et al., 2018). Manual inspection was performed to identify and remove possible isotopes left after the de-isotoping procedure. The test set peaks were matched with the common m/z vector of the cross-validation set within a search window of ±5 ppm. When multiple peaks matched, the m/z value of the most intense peak was selected. If no matched peak was found, the intensity was left equal to 0 in all the pixels. This procedure was intended to evaluate the performances of the classification model with off-sample images.

Co-localization features
A binary ROI representing the area of the image occupied by the tissue was determined by applying k-means clustering to the pre-processed intensity matrix (Inglese, Correia, et al., 2018).
After trial-and-error tests, a selection of 4 clusters was used to identify a finer segmentation of the image that could also distinguish those tissue-related areas characterized by low-intensity peaks that would cluster with the off-tissue in case of only 2 clusters. The clusters that were not localized in the corners of the image were merged to define the tissue-related ROI. This rule was determined after the visual inspection of the segmented images ( Figure 1).
Spearman's correlations between each pair of matched peaks were calculated within each image sample tissue-related ROI pixels. The asymptotic t approximation (Hollander and Wolfe, 1999) was applied to determine the significance of each correlation. The image correlations associated with a Benjamini-Hochberg corrected p-value larger than 0.05 were set equal to 0.
Subsequently, the elements of the upper triangular correlation matrix (corresponding to the correlations between different pairs of peaks) were vectorized and used as features for the image sample ( Figure 2).
No smoothing or other spatially-related intensity transformations were applied to the ion images before calculating the correlations.

Supervised modeling
The purpose of the supervised model was to predict the tumor type from the co-localization features of an off-sample mass spectrometry image. The procedure consisted of fitting a PLS-DA model on a training set and evaluation of the accuracy of the predicted labels of the validation set. 10-fold cross-validation, repeated 1,000 times, was used to assess the predictive performance of the model, using a varying number of PLS-DA components in a range of 2 to 10.
The accuracies were calculated between the original label vector and the merged predicted labels after a single repetition of the cross-validation procedure. The median accuracy over the 1,000 repetitions was considered as representative of the performances associated with the specific number of components. The maximum average accuracy determined the optimal number of PLS-DA components. A PLS-DA model fitted on the entire cross-validation set using the optimal number of components was then used to predict the labels of the test set, to assess the robustness of the predictive power of the model.
All the scripts were developed in the R language (R Development Core Team, 2016) and are available at https://github.com/paoloinglese/ion_colocalization.

RESULTS
The pre-processing workflow applied to the cross-validation set resulted in 64 common peaks. Among the 21 test set samples, one sample had one unmatched peak, two samples had two unmatched peaks, and one sample had three unmatched peaks. The 64 common peaks were annotated by generating the sum formula from the exact mass measurements and by searching the m/z values at the online Metlin database (Smith et al., 2005) (Table 3). An m/z error threshold 5 ppm was used to accept the putative annotation.
The resulting 2,016 Spearman's correlations were used as features for the PLS-DA models.
A number of 7 PLS-DA components corresponded to the maximum average accuracy across the 1,000 repeated 10-fold cross-validations, equal to 0.92050.0233 (mean, SD). A comparison with the accuracies obtained after permuting the classification labels confirmed that the observed performances were not due to random chance (Figure 3-A). The confusion matrix of a randomly selected repetition is shown in Figure 3-B.
The PLS-DA model trained on the whole cross-validation set was capable of predicting the external test set with an accuracy of 0.8095 (Figure 3-C). Although, as expected, the performance dropped, if compared with those of the cross-validation, this result confirmed that the model was adequately capable of distinguishing between tumor types.
Dependence of the correlation values and the size of the ROIs was evaluated. The number of pixels associated with the ROIs varied in a range of [1196,14345] pixels, with a median of 3726 pixels. The large number of samples (pixels) used to calculate the Spearman's correlations is compatible with previous studies of correlation stability (Schönbrodt and Perugini, 2013).
However, since in this case, the ion intensities represent samples drawn from possibly different populations (one per sub-tissue type), a further test was performed to assess the dependency of the observed correlations from the ROI sizes. Firstly, the ROI sizes (in number of pixels) were compared between the three tumor types (Figure 4-A). A Wilcoxon-Mann-Whitney test followed by Benjamini-Hochberg correction between the three pairs of ROI distributions rejected the null hypothesis of equal distribution with a significance level of 95% for the pairs 'breast-colorectal' and 'colorectal-ovarian', but not for the pair 'breast-ovarian'. In particular, breast and ovarian samples had smaller ROIs (median equal to 3205.5 and 3669.5 pixels, respectively) than colorectal samples (median equal to 9,798.5 pixels). Additionally, LOESS models between the number of ROI pixels and each correlation feature of all the cross-validation set samples showed that the average correlation values decreased with larger ROIs after (elbow point around 5,000 pixels) (Figure 4-B).
For these reasons, to test whether the classification performances were biased by the ROI size, the same classification procedure was applied to the only breast and ovarian cancer. The choice of those two classes was based on the fact that they had similar ROI size distributions with a median value lower than the elbow point of the LOESS trends.
The good classification performance of the 1,000 repeated 10-fold cross-validation (median accuracy equal to 0.8750, corresponding to 2 PLS-DA components) (Figure 4-C, D) suggested that the full model performance was not significantly biased by the different ROI sizes of the three tumor types. Analogous results were observed for the test set predictions (Figure 4-D).
A multiple pairwise comparison Dunn's test (Dunn, 1964) was employed to determine the significantly different ion-pair correlations between the tumor types. The ten most significant correlations, selected among those with more than a 95% significance level (after Benjamini-Hochberg correction) were further investigated for identifying possible hypotheses on the different biochemical mechanisms associated with the three tumor classes. Molecule pairs correlation significances were calculated as a function of the corrected p-values, defined as −10 log 10 .
The boxplots of the ten selected correlations are reported in Figure 5. The Cohen's d (Cohen, 1988) confirmed that the selected correlations had a large effect size in each pair of tumors comparisons ( Table 4).
The spatial localization of these ions could be investigated by plotting the relative abundance images of the selected correlations. The pairs of ions associated with the five most significant correlations, shown in Figure 6, confirmed that the selected ions were localized in tissue-related areas and did not involve off-tissue regions of the images.
A graph representing the top 10 ion correlations was determined from each pairwise tumor comparison. The median correlations across the samples of the same tumor class were used as the edges of a weighted graph, whereas the peak m/z corresponding to the ion pairs were assigned to the nodes. The community structure, determined via short random walks (Pons and Latapy, 2006) revealed groups of highly correlated ions among those associated with the selected correlation features. Following the assumption behind the presented workflow, these groups of ions can be used to generate hypotheses on the interactions that distinguish the biochemical characteristics of the different cancer types (Figure 7). For this purpose, the elemental composition corresponding to the m/z differences between the ion pairs associated with the selected correlations can provide information about the possible chemical interactions (Table 4).

DISCUSSION
Mass spectrometry imaging is gaining increased interest as a technology capable of mapping the spatial distribution of the molecules of interest in two or three-dimensional samples.
Unfortunately, since the quantitative analysis of imaging data is still challenging, and because of the difficulty of controlling the matrix effects hamper its use it as a quantitative tool for the determination of molecular abundances. Ion co-localization analysis, however, can represent a more robust alternative towards more quantitative analysis of this kind of data. Firstly, we showed by applying only a within-sample normalization, that co-localization patterns are not affected by the batch effect. This result was expected because Spearman's correlation detects the linear relationship between two measures independently of their range.
The workflow presented exploits the assumption that local molecular interactions are expressed through co-localization patterns, and that different classes of samples (different tumors in this case) are characterized by diverse metabolic pathways. This hypothesis was corroborated by the high predictive power of the supervised models. Similar performance was observed in an external test set, confirming the robustness of the models.
Additionally, the co-localization patterns can be used for data-driven hypothesis generation, suggesting possible local molecular mechanisms characterizing the samples of interest. The visualization of the most significantly correlated ions can help to determine regions of interest for further experiments, through data integration with transcriptomics or simply by inspection of the stained samples for histopathological validation, in the case of tissue sections. When represented as a graph the most significant correlations that reveal patterns and clusters of ions that could be used as a part of the hypothesis generation process.
Both positive and negative correlations were considered as feature values for representing the mass spectrometry images. Although a positive correlation may indicate the presence of one or multiple mechanisms locally involving the detected molecules, negative correlations may be seen as reflecting competitive processes like COX1/2, in which Arachidonic Acid is used as a common precursor whereas prostaglandins PG H2 and PG G2 are produced depending on the COX expression. Unfortunately, to date, there is currently not sufficient knowledge to derive such hypotheses exclusively the analysis of mass spectrometry imaging data.
Another limitation of the method is characterized by the representativeness of the studied samples. A single tissue section may be not representative of the complex heterogeneity of the tumor. For this reason, multiple sections or three-dimensional specimens may give a more robust representation of the molecular composition of the analyzed samples.
In the future, various similarity measures will be tested for the quantification of the degree of ion co-localization in cancer tissue and MSI data from different sources. Additionally, biologically-related hypotheses generated by the ion co-localization patterns will be thoroughly investigated through extensive experimental validation, involving transcriptomics and immunohistochemistry measurements. Inglese, P., Correia, G., et al. (2018)

Reference
Binary ROI (k-means)

Reference
Binary ROI (k-means)

Min. number of connected pixels 9
Aggressiveness 1