Feature extraction for proteomics imaging mass spectrometry data

Imaging mass spectrometry (IMS) has transformed proteomics by providing an avenue for collecting spatially distributed molecular data. Mass spectrometry data acquired with matrix assisted laser desorption ionization (MALDI) IMS consist of tens of thousands of spectra, measured at regular grid points across the surface of a tissue section. Unlike the more standard liquid chromatography mass spectrometry, MALDI-IMS preserves the spatial information inherent in the tissue. Motivated by the need to differentiate cell populations and tissue types in MALDI-IMS data accurately and efficiently, we propose an integrated cluster and feature extraction approach for such data. We work with the derived binary data representing presence/absence of ions, as this is the essential information in the data. Our approach takes advantage of the spatial structure of the data in a noise removal and initial dimension reduction step and applies $k$-means clustering with the cosine distance to the high-dimensional binary data. The combined smoothing-clustering yields spatially localized clusters that clearly show the correspondence with cancer and various noncancerous tissue types. Feature extraction of the high-dimensional binary data is accomplished with our difference in proportions of occurrence (DIPPS) approach which ranks the variables and selects a set of variables in a data-driven manner. We summarize the best variables in a single image that has a natural interpretation. Application of our method to data from patients with ovarian cancer shows good separation of tissue types and close agreement of our results with tissue types identified by pathologists.


Introduction
Mass spectrometry (MS) has become a versatile and powerful tool in proteomics for the analysis of complex biological systems, including the identification and quantification of proteins and peptides (Ong and Mann, 2005). Many different technologies have been developed under the collective field of proteomics mass spectrometry (Aebersold and Mann, 2003). The focus in this paper is the more recent development (see Groseclose et al. (2007); Aoki et al. (2007)) of matrix assisted laser desorption ionization (MALDI) imaging mass spectrometry (IMS), also known as MALDI imaging, and in particular an analysis of MALDI-IMS data acquired from tissue samples of patients with ovarian cancer.
Unlike the more common 2D gel electrophoresis (2D-GE) and liquid chromatography mass spectrometry (LC-MS) techniques in proteomics, MALDI-IMS preserves the spatial distribution inherent in the tissue; and the tens of thousands of spatially distributed spectra acquired from a single tissue sample in a MALDI-IMS experiment provide new challenges for statisticians and bioinformaticians as well as having the potential to lead to breakthroughs in biological research (see Casadonte and Caprioli (2011)). We propose a combined cluster and feature extraction method for such data which exhibits cancer specific variables whose protein associations can be inferred by parallel LC-MS experiments such as those of Meding et al. (2012).
Standard proteomics mass spectrometry methods such as 2D-GE and LC-MS have been described and reviewed in the biological, bioinformatics and statistics literature for some time, see Wasinger et al. (1995). We briefly explain LC-MS and important differences with MALDI-IMS in Section 2. For an overview and review of recent approaches in LC-MS, see America and Cordewener (2008). Statistical challenges of proteomics mass spectrometry data are outlined in Wu et al. (2003). The statistics and bioinformatics literature on the analysis of 2D-GE and LC-MS data is growing fast and covers a range of statistical methods. Testing and classification of such data are described in Morris et al. (2005); Morris (2012); Yu et al. (2006) and references therein. Other statistical approaches that have been proposed and applied to 2D-GE and LC-MS data include peak identification, alignment, and feature selection, see Yu et al. (2006); identification of proteins, see Yu et al. (2006) and Karpievitch et al. (2010); wavelet-based methods, see Morris and Carroll (2006), America and Cordewener (2008) and references therein; and methods from survival analysis for the detection of differentially expressed proteins, see Tekwe et al. (2012).
Contrasting these developments in the analyses of 2D-GE and LC-MS data, the newer MALDI-IMS methods which have been introduced into routine research practice have not yet attracted as much attention in the statistics literature, although they have been covered in proteomics/mass spectrometry journals -see Norris et al. (2007); Jones et al. (2012); Gessel et al. (2014); Stone et al. (2012) and references therein. Recent research from Alexandrov and his collaborators focuses on clustering and dimension reduction of MALDI-IMS data; see Alexandrov et al. (2010); Alexandrov and Kobarg (2011); Alexandrov et al. (2013) and references therein.
The potential of MALDI-IMS is described in Alexandrov and Kobarg (2011): 'IMS is one of the most promising innovative measurement techniques in biochemistry which has proven its potential in discovery of new drugs and cancer biomarkers. ... IMS was used in numerous studies leading to understanding chemical composition and biological processes. ... As for many modern biochemical techniques, in particular in proteomics, the development of computational methods for IMS is lagging behind the technological progress.' Related to our research are the papers by Deininger et al. (2008); Alexandrov et al. (2010) and Bonnel et al. (2011) who cluster their MALDI-IMS data using principal component analysis and hierarchical clustering, or a Gaussian mixture model, as they find that clustering of the raw data does not yield biologically interpretable results. Our proposal differs from their research in a number of important aspects. Unlike these authors, we derive suitably binned binary data instead of working with the raw or intensity data (which we describe in Section 3), and following Koch (2013, Chapter 6), who demonstrated the success in finding biologically meaningful tissue clusters for such binary data, apply k-means clustering to the binary data. Analysis of the binary data has the added advantage of being computationally more efficient since the binary data are much smaller and parallel computations are possible for all the methods we propose. Further, our approach combines clustering with explicit feature extraction -based on our Difference in Proportions of Occurrence (DIPPS) statistic which ranks and selects the 'best' variables in a data-driven way and integrates the results in a single heat map that allows a natural interpretation. The ability to visualize the selected features as a single easily interpretable image has not been a part of the above-mentioned papers and gives a significant advantage to our approach. This paper is organized in the following way. We briefly describe relevant background on proteomics in Section 2, and discuss advantages of MALDI-IMS for the biological research fields. We describe MALDI-IMS data in Section 3, indicate how to interpret such data, and derive the binary data that we will analyze. Section 4 covers the spatial smoothing and clustering steps of our method. The DIPPS approach (including a feature extraction step) is described in Section 5. In Section 6 we consider the application of our combined cluster and DIPPS approach to multiple tissue sections from three surgically excised ovarian cancers of different patients. For each cancer and each tissue section, we find a set of features that characterizes the cancer region. The results shown in Section 6 successfully separate differences between patients from within-patient variation, and illustrate how the methods can be implemented and extended in order to address biologically relevant questions such as the identification of potential tissue-specific protein markers, and classification of patients based on their response to treatment.

Proteomics Background
Ovarian cancers are virtually asymptomatic and as a result the vast majority of cases are detected when the disease is in an advanced stage and metastasis has occurred. For patients with late stage ovarian cancer, radical surgery and chemotherapy are often insufficient to address the disease adequately and many patients relapse after treatment. The combination of late-stage diagnosis and unsuccessful treatments makes ovarian cancer the most lethal gynecological cancer with only 30% of patients surviving five years following diagnosis (Ricciardelli and Oehler, 2009;Jemal et al., 2011). Compounding these clinical statistics is the heterogeneity of ovarian cancers, which makes typical proteomic analyses difficult (Deininger et al., 2008). The key to addressing ovarian cancer will be our understanding of the mechanisms driving cancer progression, the identification of molecular markers which can predict treatment success and most importantly, the identification of new treatment targets. As proteins are the main functional components of cells and tissues, determining protein distributions in cancer tissue represents a crucial step in understanding cellular dysfunction such as cancer in a holistic way.
Proteins are functional and structural components of cells. They are synthesized within cells as linear amino acid sequences and subsequently fold into more complex 3D structures that determine their functions and intra-cellular locations. The complete set of proteins which exist in a given cell, tissue or biological fluid, under defined conditions, is termed its proteome (Wilkins et al., 1996). Proteomes vary considerably between different cellular states. Understanding these differences allows for insight into the cellular processes involved in development and progression of diseases such as cancer. Proteomics is the discipline which characterizes proteome changes. Proteomic identification and quantitation, typically achieved through MS, requires the reduction of proteome complexity prior to analysis. This is often achieved with fractionation methods such as one-or two-dimensional gel electrophoresis (Gygi et al., 2000) or liquid chromatography (LC). LC is currently the predominately used fractionation technique in proteomics (Rogowska-Wrzesinska et al., 2013). In this approach molecules are affinity bound to a column material, or stationary phase, and eluted with a solvent, or mobile phase, the composition of which is usually changed as a function of time.
LC-MS analysis of a proteome uses either a top-down or bottom-up approach. Top-down approaches analyze intact proteins, whereas bottom-up approaches use proteolytic enzymes such as trypsin to digest proteins into smaller peptides prior to analysis. We are primarily interested in the bottom-up approach. In a bottom-up LC experiment peptides dissolved in a hydrophilic mobile phase are retained using hydrophobic interactions with the stationary phase. The peptides are then eluted by a linear increase in the percentage of hydrophobic solvent in the mobile phase over time. A peptide will be retained on the hydrophobic stationary phase until the peptide's affinity for the percentage organic solvent in the mobile phase exceeds its affinity for the stationary phase. The LC eluant is often directly coupled to an MS instrument, and the technique is then called LC-MS.
MS instruments contain an ion source, mass analyzer and detector. The ion source for LC-MS usually employs electro-spray ionization (ESI), but can instead use matrix assisted laser desorption ionization (MALDI) to produce gaseous peptide ions for analysis. The mass to charge ratio (m/z) of these ions is measured by the mass analyzer and detector to produce a mass spectrum, that is, measured intensity as a function of m/z. Importantly, the eluant flow from an LC system is constant, leading to the acquisition of a number of MS measurements over time. Each mass spectrum acquired in a LC-MS experiment is identified by its retention time. The peptides produced in a bottom-up proteolytic digest can also be fragmented to produce spectra which correspond to amino acid composition. Identification of these spectra is attempted using search algorithms, such as MASCOT (for more information see Koenig et al. (2008) and references therein), which match measured spectra to expected fragmentation spectra. It is the combination of LC-MS with bottom-up strategies that currently allows proteomic studies to identify and quantify thousands of unique proteins in a given biological sample.
Despite its popularity, LC-MS suffers from two method-specific limitations: 1. The length of an LC-MS run is usually measured in hours, precluding the rapid analysis of large sample numbers (≥ 20), and 2. Tissue samples used for LC-MS based proteomics are homogenized and solubilized, which removes all spatial information inherent in the tissue.
Tissues are usually a mix of different cell populations, the organization of which is directly related to their functions and hence the spatial information is essential. Different cell populations or tissue types are typically localized in distinct spatial regions, as can be seen in Figure 1(a). The loss of spatial information that occurs in LC-MS analysis has motivated the development of direct tissue analysis using MALDI-IMS Groseclose et al., 2008;Gustafsson et al., 2011).
To prepare a sample for MALDI-IMS a tissue block is thinly sectioned (2 -10µm thick slices) and mounted onto conductive microscopy slides. For an analysis of peptides, the tissue section is coated with a homogeneous layer of proteolytic enzyme and digested. The digest is overlaid with a low molecular weight organic matrix that co-crystallizes with the tissue-derived peptides. MALDI is used to ionize peptides directly from the tissue and to acquire a mass spectrum from each (x, y) position on a regular grid across the tissue section. A spacing of 20 -250µm between acquisition locations is typical and balances the requirements for high quality MS data, practical data size, and measurement time. In the ovarian cancer data a spacing of 100µm was used.
In LC-MS spectra are identified by their retention time. In contrast, MALDI-IMS spectra are curated by their spatial (x, y) coordinates, and the availability of spatial coordinates offers a unique perspective on tissue analysis. Typically, tissue structure is visualized by pathologists using haematoxylin and eosin, H&E, staining followed by light microscopy. Since tissue sections can be H&E stained post-acquisition of MALDI-IMS data (Deutskens et al., 2011), spatial changes in protein abundance measured by MALDI-IMS can be compared to the histology visible in the H&E stain of the same tissue (see Figure 1).
MALDI-IMS provides unprecedented capacities for pathology, and as a consequence has important implications for research into human cancers (Gorzolka and Walch, 2014). Pathology provides an in-depth understanding of the way in which cellular structure and appearance relate to cancer behavior and the cancer's effect on the patient. LC-MS based proteomics is able to establish identification and quantitation for large portions of the cancer proteome. MALDI-IMS, however, can combine these capacities in the following ways: 1. MALDI-IMS acquires spatially referenced data directly from cancerous tissue and has the potential to provide novel information in the form of tissue or treatment specific markers. This is not achievable with current pathology and could potentially allow clinicians to use MALDI-IMS data for individualizing patient treatments.
2. Cancers, and in particular ovarian cancers, are heterogeneous tissues (Deininger et al., 2008). Spatial analysis by MALDI-IMS allows this heterogeneity to be addressed explicitly (Gorzolka and Walch, 2014). Furthermore, if arrays of patient tumor cores are used, hundreds of patients can be analyzed in a single day. MALDI-IMS therefore provides the opportunity to screen cancer samples rapidly, target regions of interest and employ complementary LC-MS methods to characterize the tumor proteome of these tissues.
One of the key advantages of MALDI-IMS is its ability to analyze large cohorts of patients in a timely manner, for example by using a tissue micro-array or TMA (Groseclose et al., 2008;Steurer et al., 2013). To exploit this ability in future, it will be crucial to understand differences and to be able to distinguish between patients. As proof of concept, we consider data including multiple tissue types from each of three patients in this paper.

The Data and Binary Binning
In this section we describe the form of MALDI-IMS data. We also describe the steps used in transforming the initial functional data to the final 'binary binned data'. We will refer to the 'binary binned data' as 'the data' beyond this section. As discussed in Section 2, a MALDI-IMS dataset consists of mass spectra annotated by spatial (x, y) coordinates. Each mass spectrum typically contains between 5 and 200 peaks, each corresponding to a biological analyte such as a peptide or protein. A number of preprocessing steps are involved in extracting the peaks in the functional data, including: smoothing, baseline reduction, and peak picking (typically using proprietary software, e.g. flexAnalysis, Bruker Daltonik). We will only be concerned with these extracted peaks from here on.
In this paper we consider a selection of datasets from three surgically excised ovarian cancers, each from a different patient -see Gustafsson (2012). Thin crosssections of tissue are obtained from each cancer, and MALDI-IMS data are collected. We will use the term 'dataset' to refer to the data collected from a single section of tissue. Multiple sections are taken from each cancer. We expect multiple datasets from the same cancer to be similar, and aim to separate within-patient from betweenpatient variability. Initially we consider a single typifying dataset (shown in Figure 1) and will refer to it as the 'motivating dataset' throughout. A typical dataset might have 5000−100000 spectra, the motivating dataset consists of 14059. The motivating dataset will serve to illustrate our proposed method which we describe in Sections 4 and 5. In Section 6 we consider nine datasets (including the motivating dataset), three from each cancer, and compare datasets within and between patients. Although we have access to more datasets we only show nine as these suffice for comparisons both within and between patients, while still illustrating our method concisely.
As described in Section 1, each dataset contains mass spectra collected at regularly spaced points on the surface of a thin tissue section. Within a dataset, each mass spectrum is a realization of ion intensity as a function of mass-to-charge ratio (m/z). When using MALDI-MS the charge (z) is typically one and so the mass-tocharge ratio reduces to mass. Mass peaks (local maxima) can be identified and the original mass spectra reduced to 'peak list' form. This was done using the standard flexAnalysis software package (Bruker Daltonics, Bremen, Germany). Representing the data as peak lists significantly reduces the amount of data involved, often by as much as two orders of magnitude and this can be very important due to the amount of data involved (> 10GB).
The m/z locations of the peaks are discretized by binning. We call the data after binning 'binned intensity data' as in this form the variables are measures of ion intensity for each of the m/z bins. The number of variables in the binned intensity data will vary depending on the choice of bin size which turns out not to be crucial provided the chosen bin size is within a reasonable range (0.05 − 2 m/z). We will use an intermediate bin size of 0.25 m/z, which yields 4298 variables (non-empty bins) in the motivating dataset. We repeated all analyses with a number of alternate bin sizes and consistently produced very similar results, given the bin size is chosen from the above interval. The principal effect that choice of bin size has is on the total number of variables, and the number that are removed in the dimension reducing steps of our method -smaller bin sizes will result in more variables initially, and more being removed, larger bin sizes result in fewer variables initially, and fewer being removed.
Experience with exploratory analysis including k-means clustering of the binned intensity data has shown that these data do not lead to interpretable or spatially localized clusters (Koch, 2013). This is due to the effect extraneous variables (such as matrix crystal morphology) have on the intensity measurements (Garden and Sweedler, 2000). Koch (2013, Example 6.12) proposes to transform the intensity data to binary data. There exist other approaches to addressing these effects (such as the normalization approaches of Deininger et al. (2011)), but the approach of Koch (2013) bypasses the effect of the extraneous variables on peak intensities entirely. Dealing with the binary data has the added advantage of being extremely computationally efficient. The transformation to binary data is natural, as it tells us if a peak is present (takes the value one), or is absent (takes the value zero). Representing the variables in this binary form can be considered a data cleaning step, as it removes superfluous sources of noise. Subsequent analyses will refer only to the binary data, and we will refer to these simply as 'the data' from here on.

Clustering Spatially Smoothed Data
We propose a two-step method for separating spectra that belong to different tissue types. First a smoothing step, which serves dual purposes: to act as a data cleaning method, and to reduce the dimension. Second a clustering step, to separate the spectra into groups that correspond to tissue types.

Smoothing Step
The spatial information available in MALDI-IMS data can be used to clean the data and remove variables that are spatially dispersed. We incorporate this spatial information through a spatial smooth.
Let X be a d × n data matrix; the rows of X correspond to variables (m/z bins) and are denoted x i• , the columns of X correspond to spectra and are denoted x •j , and the entries of X are denoted x ij . Let 0 ≤ τ < 1 2 be a smoothing parameter and δ ≥ 0 a distance cutoff.
We iteratively update the values of X. Let X (k) denote the updated matrix at the kth iteration. Similarly let x ij denote the rows, columns, and values of X (k) respectively. At the kth iteration, the proportion of spectra, T (k) ij , in a spatial δ-neighborhood of the jth spectrum x (k) •j whose values at the ith variable x This proportion determines if the value x (k) ij should be smoothed. In (1), d j denotes a 1 × n indicator vector with entry 1 if the corresponding indexed spectrum is in a δ-neighborhood of x j• and zero otherwise. If this proportion T (k) ij is less than τ , we update the value x (k) ij as in (2). We generate the smoothed data by iteratively calculating the entries of X (k) for k = 1, 2, . . . , starting with X (0) = X. The T (k) ij term appearing in (2) is defined in (1). We stop when convergence is reached (k = k * = min{k : X (k) = X (k−1) }). The spatially smoothed data are X (k * ) .
Remarks on the smoothing process: 1. We choose δ = √ 2 which results in a range 1 Moore neighborhood. See Gray (2003). This neighborhood is used in the cellular automata literature including Conway (1970). It is worth noting that acquiring the range 1 Moore neighborhood by using the Euclidean distance and δ = √ 2 on a regular grid is equivalent to using the Tchebychev distance, and δ = 1 (without the regular grid assumption).
2. The smoothing parameter τ defines the proportion of neighboring spectra needed to agree in order for the value for an observation to remain unchanged at any given step. Small values of τ smooth less (τ = 0 leaves the data unmodified), while larger values smooth more. The limit τ → 1 2 results in maximum smoothing, and is equivalent to the intuitive median smooth. The median smooth tends to yield over-smoothed data, and often fails to converge. We chose an intermediate smoothing parameter, τ = 1 4 , for these analyses. 3. Alternative smoothing options include kernel methods (Wand and Jones, 1995) which apply to continuous data. These methods produce continuous values when applied to binary data, for which there is no clear interpretation. Our method produces binary smoothed data -maintaining the interpretability of the binary values.
4. At each smoothing iteration k, variables are smoothed independently, and within each variable all observations are smoothed simultaneously at each step. This means that it is possible to parallelize the smoothing algorithm, making efficient use of computational resources.
Variables that exhibit occurrence of peaks in a small number of spatially delocalized spectra or in almost all spectra constitute a large proportion of all variables. These variables tend not to be relevant as they are usually internal calibrants , errors, contaminants, or features that are too small (of the order 100µm) to be of interest. These variables have zero variance after the smoothing step. Following the smoothing, these zero-variance variables are removed -reducing the dimension of the data. Any such variables are not of interest as they are internal calibrants , errors, contaminants, or features that are too small (of the order 100µm) to be of interest. The motivating dataset has 4298 variables before the smoothing step. After the smoothing step 818 of these variables have non-zero variance and the remainder are discarded.

Clustering Step
The second step in our approach concerns clustering of the spatially smoothed data. We use k-means clustering with the cosine distance. Based on the information available from the histology, there are three broad tissue types present which could be labeled as fatty, connective and cancer tissue respectively. Further, there are spectra that were acquired off-tissue. Thus we perform k-means clustering with k = 4. We choose initial cluster centers as this results in the algorithm converging to a deterministic solution and allows results to be reproduced.
Even after the dimension-reducing smoothing step, the data are typically still high-dimensional. In the motivating dataset 818 variables are retained after discarding zero-variance variables at the end of the smoothing step. A natural measure of distance between high-dimensional vectors is the cosine distance -rather than the Euclidean distance, the common default in many k-means algorithms. Because we are dealing with binary data, we also repeated our analyses with the Hamming distance, which has an intuitive interpretation for binary data. The Hamming distance produces very similar results for the smoothed data, but shows slightly less stability for the unsmoothed data and so we prefer to use the cosine distance and only show the corresponding results.
For data with associated spatial information (such as MALDI-IMS data) it is natural to display the cluster membership in the form of cluster maps or cluster images: colored pixels at the (x, y) coordinates of the spectra which show the cluster membership of each spectrum using different colors to identify clusters. The H&E stained tissue cross-section and result of 4-means cluster analysis (by cosine distance) of the motivating dataset are shown in Figure 1(a) and Figure 1(b) respectively. The cluster membership in Figure 1(b) corresponds well with tissue types as determined by the histology in Figure 1(a): cyan corresponds to fatty tissue, green to cancerassociated connective tissue and orange to cancer. The cancer cluster accurately assigned the location of four cancer nodules (indicated by arrows in Figure 1(a)). The fourth cluster in pink corresponds well with off-tissue spectra. The correspondence between cluster results and histology demonstrates that the spatial smooth and cluster analysis isolate key molecular information which allows differentiation of tissue types by their mass spectra.

The Difference in Proportions of Occurrence (DIPPS) Approach
In Section 4 good agreement was demonstrated between tissue types visible in the histology of Figure 1(a) and clusters from the 4-means clustering of Figure 1(b). From the biological perspective it is of great interest to be able to quantify the differences between these tissue types in an easily interpretable way. At a mathematical level, a characterization of the differences between the tissue types translates into an identification of variables that discriminate between them, or more generally between subsets of data. We explain our approach, the DIPPS approach, for an arbitrary subset of observations with binary variables. We then apply these ideas to the motivating dataset using spectra from the cancer cluster as the subset of interest, since distinguishing these spectra from the non-cancer spectra is particularly interesting. In this section we describe the DIPPS approach which consists of 1. a ranking of the variables for their ability to characterize a subset of data; 2. a data-driven selection of the 'best' variables; and 3. a DIPPS map which combines the best variables in a single image that has a natural interpretation.

The DIPPS Approach for a General Subset of Observations
We define the DIPPS statistic for a fixed subset of the data in (3), and show how this new statistic introduces a natural ranking of the variables based on their ability to characterize the subset of interest. Further, this statistic leads to a data-driven selection of the number of ranked variables that characterize the subset. The extraction of the most highly ranked and hence most relevant variables can be interpreted as a dimension reduction step, as it effectively restricts the variables to the extracted ones. Further, we propose a way of displaying graphically the information obtained from the statistic in an easily interpretable summary image. These maps make this technique useful in exploratory analyses as combining the selected variables into a single interpretable image allows for broad conclusions to be quickly and easily drawn. When the amount of data becomes large the approach commonly used in proteomics, namely of considering each variable individually, is of limited use due to the sheer amount of information that needs to be considered simultaneously, and the ability to produce an image that summarizes all selected variables becomes particularly useful. Let S be a subset of observations of the data, i.e. a subset of the columns of X. We let p(S) denote the mean of the observations in S, and similarly let p(S c ) denote the mean of the observations in its compliment S c . Entries of the vector p(S) are the proportions of observations in S that take the value one for each variable. The interpretation of binary data as the occurrence of an event -here existence of peaks -allows the mean p(S) to be interpreted as the proportions of occurrence (for each variable) in S. These entries can be considered as the 'sensitivity' or true positive rate of each variable for S. Similarly, entries of the vector 1 d×1 − p(S c ) are the proportions of observations in S c that take the value zero for each variable and can be considered the 'specificity' or true negative rate of each variable for S. In order to characterize S, both sensitivity and specificity should be high. We sum these measures of sensitivity and specificity, and subtract one to give a range of [−1, 1] and then define the vector of DIPPS statistics for S as (3) For convenience of notation we omit the dependence on S, and write d(S) = d. We will similarly omit the S dependence for t a , n a , c, and a * below, as we are treating S as fixed. We use the entries of the vector d to rank the variables: the entry with the greatest value corresponds to the individual variable that characterizes S best. Next we determine the set of variables that collectively best characterize S. We do this by finding a cutoff value in a data-driven way as follows. For a > 0 let t a be the d × 1 vector that takes the value one if the corresponding element of d is greater than a, and is zero otherwise, and let n a be the number of ones in the vector t a . Let c be the centroid of S, calculated using the same distance, D, as in Section 4.2 -so in our case the cosine distance. Our suggested cutoff a * is where D (c, t a ) denotes the distance between c and t a . The n a * variables that best characterize S are those whose entries are one in t a * . The sum of these variables for each spectrum can then be visualized spatially in a DIPPS map.

The DIPPS Approach for the Cancer Cluster
Now we consider the DIPPS approach applied to the motivating dataset. We are particularly interested in characterizing the cancer spectra in these data, and we thus choose S to be the set of spectra belonging to the cancer cluster shown in Figure 1(b). The cutoff of (4) results in n a * = 45 variables at a * = 0.139. These variables characterize the cancer tissue, and thus warrant closer inspection. Followup analysis can be done to identify the peptides that these variables correspond to, and to draw inference as to their parent proteins. To illustrate the biological significance of these variables, we compared their m/z values to the previously published results of Gustafsson (2012). Two of the highly ranked variables (listed in Table 1, along with inferred parent proteins) were previously identified as highly expressed in cancer by manual assessment of variable distribution in datasets from this patient. The identifications were achieved through both mass matching to LC-MS/MS as well as in situ MS/MS (Gustafsson, 2012). The identity and histological distribution of these variables were successfully validated using immuno-histochemistry (unpublished results). This illustrates that the DIPPS approach can find features of known importance. Crucially, the DIPPS approach produces a list of characterizing variables more rapidly and comprehensively than a manual inspection of individual variables. An important aspect of our proposed method is the visualization of the selected variables in a single summary image which is easy to interpret, and which we call a DIPPS map. The ability to summarize results in this way is particularly attractive as considering so many variables individually can be time consuming and fail to provide an 'big-picture' perspective. We construct the DIPPS map pointwise at each (x, y) location by counting the number of extracted variables that exhibit occurrence for the spectrum at (x, y), and by visualizing this count in heat colours: cold indicating a count of zero, hot indicating that all exhibit occurrence. Figure 2(b) shows the DIPPS map for the cancer cluster.
A DIPPS map such as that shown in Figure 2(b), and the variables that define it have intuitive interpretations that the results of cluster analysis do not. This method allows variables crucial to cluster/tissue differentiation to be isolated and visualized. The visualization of such DIPPS maps becomes increasingly important in exploratory analyses when the number of patients and datasets increases as it becomes quickly infeasible to consider each of the selected variables separately. The DIPPS map highlights gradations which reveal finer detail and provides a more representative visualization of the data that is not possible in cluster maps. We discuss this point further in Section 6. Our DIPPS approach also acts as a variable reduction technique by further reducing the data from 818 to 45 variables. More importantly it successfully separates tissue type specific variables, addressing the heterogeneity of the tissue. This facilitates the comparison of datasets, which we discuss in Section 6.
6 Use of DIPPS and the Jaccard Distance in Differential Comparison of Multiple Datasets In this section we consider nine datasets, including the motivating dataset. Of these datasets three are from each of three different patients which we will refer to as patients A, B, and C, respectively. We will refer to the three datasets from patient A as A1, A2, and A3, and similarly for the datasets from patients B and C. Dataset A2 is the motivating dataset. Each of the three datasets arising from the same patient are acquired from thin (6µm) tissue sections of a single surgically excised tissue. Because of this experimental setup, we expect the cluster and DIPPS maps of datasets from the same cancer to be similar in terms of the location of the cancer clusters and the selected variables. The DIPPS map intuitively combines the variables that characterize cancer clusters, making the interpretation of large patient cohorts much easier.
The results of our analyses of the nine datasets are displayed in Figures ??, ??, and ?? for patients A, B, and C, respectively. Each figure corresponds to one cancer and shows the three datasets in rows. Each row consists of an H&E stain on the left, a cluster map in the center, and a DIPPS map (for the cancer cluster) on the right. The second row of Figure 3 repeats the results for the motivating dataset shown in Figures 1 and 2.
In all nine datasets, the clustering results correspond well with the underlying tissue morphology. Minor exceptions to this relate to dataset specific misassignment of tissue types. These exceptions do however highlight the ability of the DIPPS maps to find and extract information in the data that is not available in the cluster maps, as the DIPPS maps overcome these difficulties despite slight misalignment between cluster membership and tissue morphology. For an example of this robustness property of the DIPPS maps, consider Figure 4 -although the cluster maps show a noticeable difference in cluster shapes across the datasets from patient B, the DIPPS maps are comparatively consistent.
The DIPPS maps are also more representative of the data as they allow for gradual changes and fine detail to be seen that cannot be represented by the 'hard' boundaries in a cluster map. For example consider Figure 5, particularly C2 and C3. The DIPPS maps for these two datasets are much more representative of the cancer regions (compared with the cluster membership), as these maps indicate that the connective tissue shares variables with the cancer tissue. This is to be expected, given that the connective tissue is influenced by the presence of the cancer. However, the cancer tissue contains a higher proportion of those key cancer-specific variables, leading to a gradual change in color in the DIPPS map corresponding to fewer of these key variables occurring further away from the cancer. These gradual changes are more intuitive and informative than the cluster membership which has 'hard' boundaries and hence does not represent the data as accurately.
In addition to producing the DIPPS maps, the DIPPS approach yields the actual variables which are representative of a chosen subset. This property can be used to compare tissue types across several datasets. We calculate the DIPPS statistic and select characterizing variables as described in Section 5 for the cancer cluster in each of the nine datasets. To determine how similar these characterizing variables are, we could compare the list of variables but as the number of datasets increases, pairwise comparisons become burdensome. We use the Jaccard distance, see Jaccard (1901), for capturing similarity between subsets of observations. Other criteria could be used to calculate similarity indices, such as the χ 2 -based distances used in Koch (2013, Example 8.11) but for comparing subset overlap the Jaccard distance is a natural choice. For a pair of sets S i , S j , the Jaccard distance is We calculate the Jaccard distance for pairs of sets of characterizing variables for the cancer cluster of all nine datasets. Figure 6 shows these distances. The 'colors' range from black for a Jaccard distance of zero (indicating equality of sets) to white for a Jaccard distance of one which indicates disjoint sets. The datasets are displayed from left to right (or top to bottom) starting with those of patient A, followed by those of patient B, and finally the datasets of patient C. The block diagonal of Figure 6 is notably dark -this illustrates that differences within patients are smaller than differences between patients. This difference allows variables that are patient-specific in their characterization of the cancer to be separated from variables that are cancer specific across multiple patients. There are a number of variables that characterize cancer in all 9 of these datasets -including the m/z bin centered at 1628.75, mentioned in Section 5. The identity of this variable (listed in Table 1) has already been validated using two independent methods (in situ MS/MS and immuno-histochemistry), indicating that the protein from which this peptide derives is highly expressed in the cancer tissue of these patients. This protein could therefore be investigated further as a marker for ovarian cancer metastasis in a larger patient cohort.
On the other hand, the variables that are patient-specific in their characterization of cancer could be further investigated for their ability to classify cancers according to clinical or diagnostic criteria such as response to treatment. There appears to be a similarity between patients A and C not shared by B, as can be seen in Figure 6. This discrepancy could be a consequence of the relatively larger number of variables that characterize the cancer of patient B.' The discrepancy in variable composition between patients A/C and B can be explained by the large amount of necrotic tissue in the sections from patient B. Mass spectra from necrotic tissue are expected to be vastly different from those of other cancer tissue, which would indicate that the patient-specific variables in patient B are likely to be markers for necrotic tissue.
Other tissue types could be considered (fatty and connective tissue), and a figure similar to Figure 6 could be produced for each of them. Such figures do not show a notably darker block-diagonal in the way that Figure 6 does. This tells us that the main differences between patients are in the patients' cancer, rather than in their other tissues, and reenforces how crucial it is to address the heterogeneity of these data by separating tissue types before conducting comparisons between patients.

Conclusion
This paper proposes an integrated approach to clustering and feature extraction for functional and high-dimensional data. The research is motivated by and applied to proteomics MALDI-IMS data obtained from tissue samples and is part of a comprehensive investigation into ovarian cancer. However, our approach is generic and applies to functional or high-dimensional data with an underlying spatial structure that belong to two or more groups.
The combined cluster and feature extraction approach is based on our difference in proportions (DIPPS) statistic and includes novel visualizations which enhance the cluster maps. For the MALDI-IMS cancer data, these maps have a natural interpretation in terms of the variables that characterize cancer tissue. In other applications, the DIPPS maps overcome the limitations of single-variable visualizations and cluster maps.
A study of datasets from multiple patients shows that the DIPPS approach is able to separate variables that are patient-specific in their characterization of the cancer from variables that are cancer specific across multiple patients. This suggests that use of MALDI-IMS could lead to a differentiation of cancer types or responses of patients to potential treatments in future research. We only observed such a separation of variables for cancer tissue and not for other tissue types, which reinforced the need to address the heterogeneity of tissue data in a cluster analysis prior to feature extraction.
In proteomics, the ability to automate feature extraction and to present these features as DIPPS maps provides an opportunity for holistic appraisal of MALDI-IMS data. This is crucial due to the size and number of such datasets and their high-dimensional nature. By isolating variables important to specific tissue types and reporting similarities across patients, it will be easier to identify variables for further validation as tissue markers and to build models that could assign tissue types in new datasets by classification. Isolating variables important to specific patients and reporting differences between patients will make it easier to identify variables for further validation as predictors for addressing clinical questions such as predicting chemotherapy response and patient survival.
Yu, W., B. Wu, T. Huang, X. Li, K. Williams, and H. Zhao (2006). Statistical methods in proteomics. In Springer Handbook of Engineering Statistics, pp. 623-638. Springer. Figure 1: Left, (a) shows the H&E stained tissue section, arrows indicate the four cancer nodules. Right, (b) shows pixels plotted in their relative spatial locations, color identifies the cluster membership produced by 4-means clustering of the binary spatially smoothed data.  , and DIPPS maps for the cancer cluster (right column) for the three datasets A1, A2, and A3 (each represented in a row). Note that n a * = 43, 45, 61 variables are visualized in the DIPPS map for the datasets A1, A2, and A3, respectively.
Figure 5: H&E stains (left column), cluster maps (centre column), and DIPPS maps for the cancer cluster (right column) for the three datasets C1, C2, and C3 (each represented in a row). Note that n a * = 75, 39, 53 variables are visualized in the DIPPS map for the datasets C1, C2, and C3, respectively. Figure 6: The 9 × 9 grid shows pairwise Jaccard distances between the sets of variables characterizing the cancer clusters. Rows (and columns) correspond to the 9 datasets. The colour of each pixel indicates the value of the Jaccard distance -white corresponding to a value of one, black to zero.