Contents

1 Introduction

The goal of unsupervised analysis of mass spectrometry (MS) imaging experiments is to discover regions in the data with distinct chemical profiles, and to select the m/z-values that uniquely distinguish these different regions from each other.

Algorithmically, this means clustering the data. In imaging experiments, the resulting cluster configurations are called spatial segmentations, and the clusters are called segments.

In this vignette, we present an example segmentation workflow using Cardinal.

We begin by loading the package:

library(Cardinal)

2 Segmentation of a pig fetus wholy body cross section

This example uses the PIGII_206 dataset: a cross section of a pig fetus captured using a Thermo LTQ instrument using desorption electrospray ionization (DESI).

First, we load the dataset from the CardinalWorkflows package. The data is stored in an older format, so we need to coerce it to an MSImagingExperiment.

data(pig206, package="CardinalWorkflows")
pig206 <- as(pig206, "MSImagingExperiment")

The dataset contains 4,959 spectra with 10,200 m/z-values.

pig206
## An object of class 'MSContinuousImagingExperiment'
##   <10200 feature, 4959 pixel> imaging dataset
##     imageData(1): intensity
##     featureData(0):
##     pixelData(0):
##     run(1): PIGII_206
##     raster dimensions: 111 x 66
##     coord(2): x = 10..120, y = 1..66
##     mass range:  150.0833 to 1000.0000 
##     centroided: FALSE
Optical image of the pig fetus section

Optical image of the pig fetus section

In the optical image shown above, the brain (left), heart (center), and liver (large dark region) are clearly visible.

image(pig206, mz=885.5, plusminus=0.25)

The dataset has been cropped to remove the background slide pixels, leaving only the tissue section itself for analysis.

2.1 Pre-processing

For statistical analysis, it is useful to reduce the dataset to include only the peaks.

We calculate the mean spectrum using summarize().

pig206_mean <- summarize(pig206, .stat="mean")
plot(pig206_mean)

In order to make the mass spectra comparable between different pixels, it is necessary to normalize the data. We will use TIC normalization.

Let’s calculate the TIC to see how it currently varies across the dataset in the raw, unprocessed specra.

pig206_tic <- summarize(pig206, .stat=c(tic="sum"), .by="pixel")
image(pig206_tic)

To process the dataset, we will first perform peak picking on the mean spectrum to create a set of reference peaks. We will then bin the peaks in the entire dataset to this reference.

Note that peak picking on the mean spectrum is the fastest option, but may miss low-intensity peaks or peaks that only occur in a small part of the dataset. If we wanted to be more thorough, we could use a similar procedure to perform peak picking on the entire dataset (or on a random sample of many spectra) to create the set of reference peaks.

pig206_ref <- pig206_mean %>%
  peakPick(SNR=3) %>%
  peakAlign(ref="mean",
            tolerance=0.5,
            units="mz") %>%
  peakFilter() %>%
  process()
## nrows changed from 10200 to 812

Now we bin the rest of the dataset to the reference peaks.

pig206_peaks <- pig206 %>%
  normalize(method="tic") %>%
  peakBin(ref=mz(pig206_ref),
          tolerance=0.5,
          units="mz") %>%
  process()

pig206_peaks
## An object of class 'MSContinuousImagingExperiment'
##   <110 feature, 4959 pixel> imaging dataset
##     imageData(1): intensity
##     featureData(0):
##     pixelData(0):
##     processing complete(2): normalize peakBin
##     processing pending(0):
##     run(1): PIGII_206
##     raster dimensions: 111 x 66
##     coord(2): x = 10..120, y = 1..66
##     mass range: 157.3648 to 889.5534 
##     centroided: TRUE

This produces a centroided dataset with 110 peaks.

2.2 Visualization

Before proceeding with the statistical analysis, we’ll first perform some and exploratory visual analysis of the dataset.

2.2.1 Ion images

Below, we plot several hand-selected peaks corresponding to major organs.

m/z 187 appears highly abundant in the heart.

image(pig206_peaks, mz=187)

m/z 840 appears highly abundant in the brain and spinal cord.

image(pig206_peaks, mz=840)

m/z 537 appears highly abundant in the liver.

image(pig206_peaks, mz=537)

Rather than manually going the full dataset and hand-selecting peaks, the goal of our statistical analysis will be to automatically select the peaks that distinguish such regions (e.g., the major organs).

2.2.2 Principal components analysis (PCA)

Principal component analysis (PCA) is a popular method for exploring a dataset. PCA is available in Cardinal through the PCA() method.

Below, we calculate the first 3 principal components.

pig206_pca <- PCA(pig206_peaks, ncomp=3)

Next, we overlay the first 3 principal components.

The overlay requires some contrast enhancement to see the structures clearly. In addition, the range of the PC scores are normalized to the same range (between 0 and 1).

image(pig206_pca, contrast.enhance="histogram", normalize.image="linear")

We can plot the loadings for the principal components as well.

plot(pig206_pca, lwd=2)

PCA can sometimes be useful for exploring a dataset. For example, here, we can see that PC3 appears to distinguish the liver, but also includes several other structures. This makes it difficult to fully utilize PCA for analysis.

2.3 Segmentation with spatial shrunken centroids (SSC)

To segment the dataset and automatically select peaks that distinguish each region, we will use the spatialShrunkenCentroids() method provided by Cardinal.

Important parameters to this method include:

  • method The type of spatial weights to use:

    • “gaussian” weights use a simple Gaussian smoothing kernel

    • “adaptive” weights use an adaptive kernel that sometimes preserve edges better

  • r The neighborhood smoothing radius; this should be selected based on the size and granularity of the spatial regions in your dataset

  • s The shrinkage or sparsity parameter; the higher this number, the fewer peaks will be used to determine the final segmentation.

  • k The maximum number of segments to try; empty segments are dropped, so the resulting segmentation may use fewer than this number.

It is typically best to set k relatively high and let the algorithm drop empty segments. You typically also want to try a wide range of sparsity with the s parameter.

set.seed(1)
pig206_ssc <- spatialShrunkenCentroids(pig206_peaks, method="adaptive",
                                       r=2, s=c(0,5,10,15,20,25), k=10)
summary(pig206_ssc)
## Spatially-aware nearest shrunken centroids:
##  
##  Segmentation / clustering 
##  Method = adaptive 
##  Distance = chebyshev
##  
##   Radius (r) Init (k) Shrinkage (s) Classes Features/Class
## 1          2       10             0      10         110.00
## 2          2       10             5      10          67.90
## 3          2       10            10       9          40.56
## 4          2       10            15       6          34.67
## 5          2       10            20       6          22.50
## 6          2       10            25       5          15.60

As shown in the summary, the resulting number of segments typically decreases as s increases. This is because fewer peaks are used to determine the segmentation.

First, non-informative peaks are removed, but as s increases meaningful peaks may be removed as well. The most interesting and useful segmentations tend to be the ones with the highest value of s just before the resulting number of segments decreases too much.

2.3.1 Plotting the segmentation

Let’s plot the resulting segmentations for s = 10, 15, 20, 25.

image(pig206_ssc, model=list(s=c(10,15,20,25)))

It is useful to see how the segmentation changes as fewer peaks are used and the number of segments decreases. Noisy, less-meaningful segments tend to be removed first, so we want to explore the segmentation with the highest value of s that still captures most of the regions we would expect to see.

image(pig206_ssc, model=list(s=20))

Here, we can see the heart, brain, and liver distinguished as segments 1, 5, and 6.

2.3.2 Plotting the (shrunken) mean spectra

Plotting the shrunken centroids is analogous to plotting the mean spectrum of each segment.

plot(pig206_ssc, model=list(s=20), lwd=2)

Let’s break out the centroids for the heart, brain, and liver segments.

cols <- discrete.colors(6)
setup.layout(c(3,1))
plot(pig206_ssc, model=list(s=20), column=1, col=cols[1], lwd=2, layout=NULL)
plot(pig206_ssc, model=list(s=20), column=5, col=cols[5], lwd=2, layout=NULL)
plot(pig206_ssc, model=list(s=20), column=6, col=cols[6], lwd=2, layout=NULL)

Some differences are visible, but it can be difficult to tell exactly which peaks are changing between different segments based on the mean spectra alone.

2.3.3 Plotting and interpretting t-statistics of the m/z values

Plotting the t-statistics tells us exactly the relationship between each segment’s centroid and the global mean spectrum. The t-statistics are the difference between a segment’s centroid and the global mean, divided by a standard error.

Positive t-statistics indicate that peak is systematically higher in that segment as compared to the global mean spectrum.

Negative t-statistics indicate that peak is systematically lower in that segment as compared to the global mean spectrum.

Spatial shrunken centroids works by shrinking these t-statistics toward 0 by s, and using the new t-statistics to recompute the segment centroids. The effect is that peaks that are not very different between a specific segment and the global mean are effectively eliminated from the segmentation.

plot(pig206_ssc, model=list(s=20), values="statistic", lwd=2)

If we break out the t-statistics for the heart, brain, and liver segments we can learn something interesting.

setup.layout(c(3,1))
plot(pig206_ssc, model=list(s=20), values="statistic",
     column=1, col=cols[1], lwd=2, layout=NULL)
plot(pig206_ssc, model=list(s=20), values="statistic",
     column=5, col=cols[5], lwd=2, layout=NULL)
plot(pig206_ssc, model=list(s=20), values="statistic",
     column=6, col=cols[6], lwd=2, layout=NULL)

Very few peaks distinguish the heart, while many more distinguish the brain and liver.

2.3.4 Retrieving the top m/z-values

Use the topFeatures() method to extract the m/z values of the peaks that most distinguish each segment, ranked by t-statistic.

Peaks associated with the heart:

topFeatures(pig206_ssc, model=list(s=20), class==1)
## Top-ranked features: 
##          mz r  k  s class   centers statistic
## 1  187.3026 2 10 20     1 315.61657 44.756915
## 2  179.3393 2 10 20     1  98.81549  4.936801
## 3  157.3648 2 10 20     1  13.56694  0.000000
## 4  159.3490 2 10 20     1  23.59401  0.000000
## 5  166.2638 2 10 20     1  34.16502  0.000000
## 6  167.2917 2 10 20     1  28.88189  0.000000
## 7  169.2804 2 10 20     1  17.83297  0.000000
## 8  171.3623 2 10 20     1  11.82048  0.000000
## 9  174.3602 2 10 20     1  19.70869  0.000000
## 10 175.3548 2 10 20     1  36.83639  0.000000

Peaks associated with the brain:

topFeatures(pig206_ssc, model=list(s=20), class==5)
## Top-ranked features: 
##          mz r  k  s class   centers statistic
## 1  885.5497 2 10 20     5 136.80283  41.60596
## 2  834.3988 2 10 20     5  62.54057  35.96851
## 3  840.4124 2 10 20     5  64.81401  34.20209
## 4  812.4232 2 10 20     5  75.69838  28.53206
## 5  305.5080 2 10 20     5  67.06297  28.50655
## 6  886.5588 2 10 20     5  65.70995  28.11709
## 7  838.4014 2 10 20     5  40.02901  27.53208
## 8  810.4166 2 10 20     5  55.01009  26.51290
## 9  887.5696 2 10 20     5 136.55897  26.11239
## 10 303.4838 2 10 20     5  80.21382  22.57846

Peaks associated with the liver:

topFeatures(pig206_ssc, model=list(s=20), class==6)
## Top-ranked features: 
##          mz r  k  s class   centers statistic
## 1  536.9240 2 10 20     6 174.22883  96.97570
## 2  562.9183 2 10 20     6 161.47592  95.65799
## 3  534.9808 2 10 20     6 119.05397  83.20395
## 4  887.5696 2 10 20     6 199.38260  48.60206
## 5  885.5497 2 10 20     6 164.42736  47.36604
## 6  281.5297 2 10 20     6 301.31716  42.23976
## 7  508.9778 2 10 20     6  55.88272  39.40409
## 8  253.5862 2 10 20     6  85.87183  37.40256
## 9  327.2898 2 10 20     6 101.12881  37.20571
## 10 888.5543 2 10 20     6  90.32627  33.71810

The top m/z values for each segment match up well with the hand-selected peaks.

3 Segmentation of a cardinal painting

It can be difficult to evaluate unsupervised methods (like segmentation) on data where we do not know the ground truth.

In this section, we use an MS image of a painting, where we know the ground truth.

data(cardinal, package="CardinalWorkflows")
cardinal <- as(cardinal, "MSImagingExperiment")
Cardinal painting

Cardinal painting

In this experiment, DESI spectra were collected from an oil painting of a cardinal.

cardinal
## An object of class 'MSContinuousImagingExperiment'
##   <10800 feature, 12600 pixel> imaging dataset
##     imageData(1): intensity
##     featureData(0):
##     pixelData(0):
##     run(1): Bierbaum_demo_
##     raster dimensions: 120 x 105
##     coord(2): x = 1..120, y = 1..105
##     mass range:  100.0833 to 1000.0000 
##     centroided: FALSE

The dataset includes 12,600 spectra with 10,800 m/z values.

3.1 Pre-processing

First, we will pre-process the dataset as before, by applying peak picking to the mean spectrum.

cardinal_mean <- summarize(cardinal, .stat="mean")
cardinal_ref <- cardinal_mean %>%
  peakPick(SNR=3) %>%
  peakAlign(ref="mean",
            tolerance=0.5,
            units="mz") %>%
  peakFilter() %>%
  process()
## nrows changed from 10800 to 1183
cardinal_peaks <- cardinal %>%
  normalize(method="tic") %>%
  peakBin(ref=mz(cardinal_ref),
          tolerance=0.5,
          units="mz") %>%
  process()

cardinal_peaks
## An object of class 'MSContinuousImagingExperiment'
##   <106 feature, 12600 pixel> imaging dataset
##     imageData(1): intensity
##     featureData(0):
##     pixelData(0):
##     processing complete(2): normalize peakBin
##     processing pending(0):
##     run(1): Bierbaum_demo_
##     raster dimensions: 120 x 105
##     coord(2): x = 1..120, y = 1..105
##     mass range: 101.0496 to 650.1034 
##     centroided: TRUE

This results in a centroided dataset with 106 peaks.

3.2 Segmetation with SSC

Now we use spatial shrunken centroids to segment the dataset.

set.seed(1)
cardinal_ssc <- spatialShrunkenCentroids(cardinal_peaks, method="adaptive",
                                       r=2, s=c(10,20,30,40), k=10)
summary(cardinal_ssc)
## Spatially-aware nearest shrunken centroids:
##  
##  Segmentation / clustering 
##  Method = adaptive 
##  Distance = chebyshev
##  
##   Radius (r) Init (k) Shrinkage (s) Classes Features/Class
## 1          2       10            10      10          36.20
## 2          2       10            20       8          21.75
## 3          2       10            30       7          14.14
## 4          2       10            40       7          10.00
image(cardinal_ssc)

We can see that with s = 10 and s = 20, two segmments are capturing an unwanted background gradient. At s = 30, this background gradient is eliminated.

Now we can use the segmentation to re-construct the original painting.

image(cardinal_ssc, model=list(s=40),
      col=c("1"=NA, "2"="gray", "3"="black", "4"="firebrick",
            "5"="brown", "6"="darkred", "7"="red"))

Let’s find the m/z values associated with the cardinal’s body.

topFeatures(cardinal_ssc, model=list(s=40), class==7)
## Top-ranked features: 
##          mz r  k  s class   centers statistic
## 1  207.0370 2 10 40     7 853.41434 206.41178
## 2  277.1958 2 10 40     7 543.28804 189.91259
## 3  291.1863 2 10 40     7 413.34428 170.52895
## 4  418.6775 2 10 40     7 411.98214 102.22334
## 5  327.1961 2 10 40     7 224.95026  98.55784
## 6  305.1687 2 10 40     7 106.09898  68.77804
## 7  208.0414 2 10 40     7  60.22322  52.26126
## 8  278.2111 2 10 40     7  49.66764  42.97868
## 9  419.7390 2 10 40     7  43.24051  28.87697
## 10 292.2195 2 10 40     7  32.31241  28.09562
image(cardinal, mz=207)

And let’s find the m/z values associated with the “DESI-MS” text.

topFeatures(cardinal_ssc, model=list(s=40), class==3)
## Top-ranked features: 
##          mz r  k  s class    centers statistic
## 1  649.0976 2 10 40     3 1053.52694 138.74473
## 2  650.1034 2 10 40     3  323.45026  94.97967
## 3  101.0496 2 10 40     3   42.72138   0.00000
## 4  112.9746 2 10 40     3   51.05414   0.00000
## 5  115.0517 2 10 40     3  112.25387   0.00000
## 6  119.9959 2 10 40     3   18.94445   0.00000
## 7  121.0219 2 10 40     3   14.46792   0.00000
## 8  127.0095 2 10 40     3   30.19322   0.00000
## 9  129.0518 2 10 40     3   74.26432   0.00000
## 10 133.9885 2 10 40     3   39.39943   0.00000
image(cardinal, mz=649)

4 Session information

sessionInfo()
## R version 3.6.1 (2019-07-05)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.3 LTS
## 
## Matrix products: default
## BLAS:   /home/biocbuild/bbs-3.10-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.10-bioc/R/lib/libRlapack.so
## 
## Random number generation:
##  RNG:     L'Ecuyer-CMRG 
##  Normal:  Inversion 
##  Sample:  Rejection 
##  
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=C              
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats4    parallel  stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
## [1] CardinalWorkflows_1.18.0 Cardinal_2.4.0          
## [3] ProtGenerics_1.18.0      S4Vectors_0.24.0        
## [5] EBImage_4.28.0           BiocParallel_1.20.0     
## [7] BiocGenerics_0.32.0      BiocStyle_2.14.0        
## 
## loaded via a namespace (and not attached):
##  [1] biglm_0.9-1        locfit_1.5-9.1     tidyselect_0.2.5  
##  [4] xfun_0.10          purrr_0.3.3        lattice_0.20-38   
##  [7] htmltools_0.4.0    viridisLite_0.3.0  yaml_2.2.0        
## [10] rlang_0.4.1        pillar_1.4.2       glue_1.3.1        
## [13] DBI_1.0.0          sp_1.3-1           jpeg_0.1-8.1      
## [16] stringr_1.4.0      htmlwidgets_1.5.1  evaluate_0.14     
## [19] Biobase_2.46.0     knitr_1.25         irlba_2.3.3       
## [22] Rcpp_1.0.2         BiocManager_1.30.9 abind_1.4-5       
## [25] png_0.1-7          digest_0.6.22      stringi_1.4.3     
## [28] tiff_0.1-5         bookdown_0.14      dplyr_0.8.3       
## [31] grid_3.6.1         tools_3.6.1        bitops_1.0-6      
## [34] magrittr_1.5       RCurl_1.95-4.12    tibble_2.1.3      
## [37] crayon_1.3.4       pkgconfig_2.0.3    MASS_7.3-51.4     
## [40] Matrix_1.2-17      matter_1.12.0      assertthat_0.2.1  
## [43] rmarkdown_1.16     R6_2.4.0           fftwtools_0.9-8   
## [46] mclust_5.4.5       signal_0.7-6       nlme_3.1-141      
## [49] compiler_3.6.1