Morphometric Cell Classification for Single‐Cell MALDI‐Mass Spectrometry Imaging

Abstract The large‐scale and label‐free molecular characterization of single cells in their natural tissue habitat remains a major challenge in molecular biology. We present a method that integrates morphometric image analysis to delineate and classify individual cells with their single‐cell‐specific molecular profiles. This approach provides a new means to study spatial biological processes such as cancer field effects and the relationship between morphometric and molecular features.


MALDI-MSI data processing
Porcine colon dataset: The MALDI-timsTOF fleX data (130,479 pixels, resolving power 40 000 at m/z 800) in its SQLite format was loaded into MATLAB R2018b (MathWorks, USA). The spectrum from each pixel was normalized to its root mean square intensity value. The overall mean spectrum from all pixels was calculated and subjected to peak picking using the MATLAB command mspeaks (Bioinformatics toolbox v2018a) with a minimum absolute peak intensity threshold of 5. This commands returns detected lipid masses, their intensities and peak extensions after default wavelet de-noising.
The resulting peak list contained 1252 m/z intervals (supplementary figure 5). These intervals were used to reduce the original data by extracting the maximum intensity within the extensions of each interval from every spectrum, resulting in a 130,479 × 1252 data matrix.
Diffuse-type gastric carcinoma dataset: The raw MALDI-TOF rapifleX data (47,355 pixels, resolving power 20 000 at m/z 800) was exported as imzML format from FlexImaging (v5.0, Bruker Daltonik) and loaded into MATLAB R2018b (MathWorks). The spectrum from each pixel was normalized to its total ion count. The overall mean spectrum from all pixels was calculated and subjected to peak picking using the MATLAB command mspeaks (Bioinformatics toolbox v2018a) with a minimum absolute peak intensity threshold of 10E-4. This command returns detected lipid masses, their intensities and peak extensions after wavelet de-noising. The resulting peak list contained 334 m/z intervals (supplementary figure 12). These intervals were used to reduce the original data by extracting the maximum intensity within the extensions of each interval from every spectrum, resulting in a 47,355 × 334 data matrix.

Histological staining & slide scanning
After MSI experiments, all slides were washed in 100% methanol (Merck KGaA, Germany) to remove the matrix and subjected to conventional hematoxylin & eosin (H&E) histological staining. For this, the slides were submersed for 3 min in each of the series of ethanol (Merck KGaA) solutions: 100%, 100%, 96%, 96% and 70%, 70%. Next, the slides were washed for 3 min in water and submersed for 3 min in hematoxylin (Merck KGaA). Then, the hematoxylin solution was washed off under running water, and the slides were submersed in eosin (Merck KGaA) for 30 sec. The eosin was washed off under running water and the slides were submersed in 100% ethanol for 1 min and 100% xylene (Merck KGaA) for 5 min. Finally, the slides were mounted on glass cover slips using Entellan mounting medium (Merck KGaA) followed by drying for minimally 12 h.
After drying, each slide was scanned by a digital slide scanner (Aperio CS2, Leica, Germany) at 20× magnification. The

Morphometric image analysis & cell classification
The scanned high-resolution histological image was loaded into the QuPath software (v0.2.0.m8) where the following steps were carried out: 1. Manual definition of regions of interest (ROIs) for each cell type to be classified: a. Porcine colon dataset: 24 ROIs for glandular cells and 17 ROIs for lamina propria cells.
b. Diffuse-type gastric carcinoma (DGC) dataset: 5 ROIs for tumor cells and 4 ROIs for muscle cells.
2. Automated cell detection in the defined ROIs using the following parameters:
The intensity-based, automatic registration between the processed binary images was performed using the multimodal configuration and affine transformation. This registration resulted in a correction of the H&E image, represented as ΔHE-HE_corrected.
In the porcine colon dataset, the resulting error could be assessed by measuring the Euclidean distance between the H&E pixel positions provided by the geometrical transformation and the actual visible laser burns (Supplementary Figure 6).

Linking single cells to their molecular profiles
Once MALDI-MSI and H&E images were aligned, the annotations from QuPath were imported from the previously exported text file. The coordinates of each cell were then transformed by first subtracting the crop points xmin and ymin and subsequently applying the geometric transformation ΔHE-HE_corrected.
To be able to overlay the MSI raster with the cell circumferences, for every MSI pixel pi with coordinates pi,x and pi,y four point coordinates were created to describe its square area, done by subtracting and/or adding half the MSI pixel size (here 5 µm) to pi: [(pi,x-5, pi,y+5); (pi,x-5, pi,y-5); (pi,x+5, pi,y-5); (pi,x+5, pi,y+5)]. These MSI square coordinates were transformed to the H&E space by first applying the geometric transformation ΔHE-MSI and then by subtracting the crop points xmin and ymin.
Every cell was then checked for its areal overlap with the areas of the 15 closest MSI pixels. The cell's molecular profile was calculated using a weighted-mean across all 15 MSI pixels, where the weights correspond to the degree of overlap between the cell's area and each MSI pixel's area, representing the contribution of each MSI pixel to the cell's molecular profile.
Scheme 1: Visual representation of the weighted average strategy to compute single-cell molecular profiles from several MSI pixels.

Spatial statistics
The centroid position of each cell was calculated and used to compute a distance matrix containing the geometric separation between all single cells. This distance matrix was used in the DGC sample to compare the molecular profiles given by MALDI-MSI between muscle cells close to tumor cells and muscle cells far away from tumor cells. We chose to only to take into account and sort distances between muscle cells and tumor cells. An average proximity of a given muscle cell to its k-nearest tumor cells was calculated, where k was chosen to be 5. These average proximities were thresholded at the first and third quartiles to define muscle cells close to and far away from tumor cells, respectively.

Correlation analysis of morphometric features with molecular features
First, a principal component analysis (PCA) was performed in MATLAB R2018b using statistics toolbox on the 47,355 pixels × 334 peaks data matrix of the DGC MALDI-MSI dataset. As described above, the principal component scores of the PCA were averaged per cell based on the areal overlap of the cell to the MSI pixels. Next, the morphometric features detected for every cell were exported from QuPath as CSV file and imported to MATLAB. Finally, a Pearson Spearman correlation analysis was performed between the first five principal components of the PCA and the 41 morphometric features.