Contents

1 Introduction

For experiments in which analyzed samples come from different classes or conditions, a common goal of statistical analysis is class comparison via hypothesis testing.

Statistical testing is performed to find peaks that are differentially abundant among different classes or conditions.

Valid statistical testing requires biological replicates in order to compare between different conditions. It should never be performed with less than 3 samples per condition.

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

We begin by loading the package:

library(Cardinal)

2 Statistical testing for a renal cell carcinoma (RCC) cancer dataset

This example uses DESI spectra collected from a renal cell carcinoma (RCC) cancer dataset consisting of 8 matched pairs of human kidney tissue. Each tissue pair consists of a normal tissue sample and a cancerous tissue sample.

MH0204_33 UH0505_12 UH0710_33 UH9610_15
UH9812_03 UH9905_18 UH9911_05 UH9912_01

For the RCC cancer dataset, the goal is to find m/z features differentially abundant between normal and cancer tissue.

First, we load the dataset from the CardinalWorkflows package using exampleMSIData().

rcc <- CardinalWorkflows::exampleMSIData("rcc")

The dataset contains 16,000 spectra with 10,200 m/z-values.

rcc
## MSImagingExperiment with 10200 features and 16000 spectra 
## spectraData(1): intensity
## featureData(1): mz
## pixelData(4): x, y, run, diagnosis
## coord(2): x = 1...99, y = 1...38
## runNames(8): MH0204_33, UH0505_12, UH0710_33, ..., UH9905_18, UH9911_05, UH9912_01
## mass range:  150.08 to 1000.00 
## centroided: FALSE

We can visualize the diagnoses:

image(rcc, "diagnosis", layout=c(2,4), free="xy", col=dpal("Set1"))

As can be seen above, each matched pair of tissues belonging to the same subject are on the same slide (i.e., the same run). Note also the the cancer tissue is on the left and the normal tissue is on the right on each slide.

2.1 Pre-processing

First, let’s visualize the total ion current.

rcc <- summarizePixels(rcc, stat=c(TIC="sum"))
image(rcc, "TIC", layout=c(2,4), free="xy")

Clearly there is pixel-to-pixel variation in in the total ion current.

We will normalize the TIC, perform peak picking on a sample of 10% of all mass spectra to create a set of reference peaks, and then summarize these reference peaks in every spectrum.

rcc_peaks <- rcc |>
    normalize(method="tic") |>
    peakProcess(SNR=3, sampleSize=0.1, filterFreq=0.2,
        tolerance=0.5, units="mz")

rcc_peaks
## MSImagingExperiment with 266 features and 16000 spectra 
## spectraData(1): intensity
## featureData(3): mz, count, freq
## pixelData(5): x, y, run, diagnosis, TIC
## coord(2): x = 1...99, y = 1...38
## runNames(8): MH0204_33, UH0505_12, UH0710_33, ..., UH9905_18, UH9911_05, UH9912_01
## metadata(1): processing
## mass range: 151.2379 to 994.4166 
## centroided: TRUE

Rather than rely on the manual region-of-interest selection, we will rely on the fact that cancer tissue is on the left and the normal tissue is on the right on each slide.

x_threshold <- c(35, 23, 28, 39, 29, 28, 47, 32)

rcc_peaks$rough_diagnosis <- factor("normal", levels=c("cancer", "normal"))

for ( i in seq_len(nrun(rcc_peaks)) ) {
    irun <- run(rcc_peaks) == runNames(rcc_peaks)[i]
    j <- irun & coord(rcc_peaks)$x < x_threshold[i]
    pData(rcc_peaks)$rough_diagnosis[j] <- "cancer"
}

rcc_peaks$samples <- interaction(run(rcc_peaks), rcc_peaks$rough_diagnosis)

We can plot the split to double-check that it is correct.

rcc_peaks$cancer <- ifelse(rcc_peaks$rough_diagnosis=="cancer",
    rcc_peaks$TIC, 0)
rcc_peaks$normal <- ifelse(rcc_peaks$rough_diagnosis=="normal",
    rcc_peaks$TIC, 0)

image(rcc_peaks, c("cancer", "normal"), superpose=TRUE,
    layout=c(2,4), free="xy", col=dpal("Set1"),
    enhance="histogram", scale=TRUE)

2.1.1 Non-specific filtering to reduce data size

In order to reduce the size of the dataset further (because the computation we are working toward can be time consuming), we will perform non-specific filtering.

This means filtering our peaks based on a summary statistic unrelated to the condition. We will use the variance.

rcc_peaks <- summarizeFeatures(rcc_peaks, stat=c(Variance="var"))

plot(rcc_peaks, "Variance", xlab="m/z", ylab="Intensity")

Now we keep only the peaks above the top 80% quantile of variance among peaks.

rcc_peaks <- subsetFeatures(rcc_peaks, Variance >= quantile(Variance, 0.8))

rcc_peaks
## MSImagingExperiment with 54 features and 16000 spectra 
## spectraData(1): intensity
## featureData(4): mz, count, freq, Variance
## pixelData(9): x, y, run, ..., samples, cancer, normal
## coord(2): x = 1...99, y = 1...38
## runNames(8): MH0204_33, UH0505_12, UH0710_33, ..., UH9905_18, UH9911_05, UH9912_01
## metadata(1): processing
## mass range: 157.3255 to 887.4361 
## centroided: TRUE

This reduces the dataset to only 54 peaks.

2.1.2 Segmentation with spatial Dirichlet Gaussian mixture model (DGMM)

Spatial-DGMM performs spatial segmentation on every peak in the dataset, and estimates the means and variances of the segments.

This is useful because we can use it to automatically detect segments to compare for statistical testing (e.g., “cancer” vs “normal” tissue). However, to do this in a way that doesn’t bias downstream statistical testing, we must make sure the segmentation is performed independently for each sample. We do this by specifying the tissue samples to the groups argument.

set.seed(1)
rcc_dgmm <- spatialDGMM(rcc_peaks, r=1, k=3, groups=rcc_peaks$samples)

rcc_dgmm
## SpatialDGMM on 54 variables and 16000 observations
## names(10): class, mu, sigma, ..., weights, r, k
## coord(2): x = 1...99, y = 1...38
## runNames(8): MH0204_33, UH0505_12, UH0710_33, ..., UH9905_18, UH9911_05, UH9912_01
## modelData(): Spatial Gaussian mixture models (k=3, n=54)
## 
## Parameter estimates:
## $mu
## , , 1 
##                           1          2          3
## MH0204_33.cancer 28.4431187 17.0674768  0.7647914
## UH0505_12.cancer 15.9890383 11.1571133  1.9363772
## UH0710_33.cancer 20.7573350 13.3081602  1.9170453
## UH9610_15.cancer 25.0320739 15.9287176  4.6695475
## UH9812_03.cancer 16.7212698  8.9047147  0.7593033
## UH9905_18.cancer 15.4596126 11.8858005  6.7703981
## ...                     ...        ...        ...
## , , ... 
## 
## $sigma
## , , 1 
##                          1         2         3
## MH0204_33.cancer 4.6953330 3.5727040 1.8632391
## UH0505_12.cancer 1.5720243 1.8009069 1.8302321
## UH0710_33.cancer 3.0549486 2.6017867 1.9436316
## UH9610_15.cancer 4.0631834 3.4427290 2.2436886
## UH9812_03.cancer 3.4825164 2.1569670 1.4231585
## UH9905_18.cancer 1.7393026 0.9772716 2.2529100
## ...                    ...       ...       ...
## , , ...

We can see from the estimated parameters above that each sample was segmented independently and have their own parameter estimates.

After segmenting the images for every peak, the segments are ranked according to their mean intensities. We can then visualize their ranks.

image(rcc_dgmm, i=2, layout=c(2,4), free="xy")

We can also plot the estimated parameters for each segment for a peak.

plot(rcc_dgmm, i=2)

Each curve shows the estimated intensity distribution for each segment (of 3) and each sample (of 16). Curves of the same color represent segments of the same rank in different samples.

2.2 Class comparison with means-based testing

As introduced earlier, statistical testing is performed to find peaks differentially abundant among different groups. Since MS imaging produces many hundreds of measurements on the same sample, we can’t treat each mass spectrum as a separate observation. Rather, we need to compare entire samples rather than individual pixels.

One way to do this is to summarize each sample by calculating its mean intensity. We can then fit linear models to the means-summarized data.

2.2.1 Fitting models with means-summarized groups

In Cardinal, we can simply use meansTest() to perform means-based testing in a MS imaging experiment. We use a one-sided formula to specify the fixed effects (the diagnosis in this particular dataset). The samples must also be provided. Each sample is summarized by its mean, and then a linear model is fit to the summaries.

rcc_mtest <- meansTest(rcc_peaks, ~rough_diagnosis, samples=rcc_peaks$samples)

rcc_mtest
## MeansTest of length 54
## model: lm 
##   i       mz                       fixed   statistic      pvalue
## 1 1 157.3255 intensity ~ rough_diagnosis  0.43003571 0.511971425
## 2 2 171.3083 intensity ~ rough_diagnosis  0.96326841 0.326364775
## 3 3 179.1700 intensity ~ rough_diagnosis  0.05159943 0.820303164
## 4 4 181.1523 intensity ~ rough_diagnosis  0.27239479 0.601729834
## 5 5 201.2486 intensity ~ rough_diagnosis  0.49512497 0.481650037
## 6 6 215.2422 intensity ~ rough_diagnosis 10.22674596 0.001384187
## ... and 48 more results

2.2.2 Interpreting the results

We can use the topFeatures() method to find differentially abundant peaks. The ranked results are automatically adjusted for multiple comparisons using the false discovery rate, or FDR. We now look for tests that are statistically significant at the 0.05 level.

rcc_mtest_top <- topFeatures(rcc_mtest)

subset(rcc_mtest_top, fdr < 0.05)
## DataFrame with 0 rows and 5 columns

But we don’t find any.

2.3 Class comparison with segmentation-based testing

Means-based testing is fast and simple and can work well for homogeneous samples. However, it doesn’t use the spatial structure of each peak, so it doesn’t take advantage of MS imaging, and may result low statistical power (i.e., failing to detect differences that actually exist).

Rather than simply average the intensities, we can summarize each sample by segmenting it with spatial-DGMM, and comparing the resulting segments. This gives us a bias-free way to keep the spatial heterogeneous information.

2.3.1 Fitting models with spatial-DGMM-summarized groups

First, we must segment the data with spatialDGMM(), while making sure that each sample is segmented independently (as specified by groups). We’ve already done this. Now we use meansTest() on the resulting SpatialDGMM segmentation to fit the models.

The segments within each sample are ranked by their mean intensities, and the segments with the top-ranked means are used for comparison.

rcc_stest <- meansTest(rcc_dgmm, ~rough_diagnosis)

rcc_stest
## MeansTest of length 54
## model: lm 
##   i       mz                       fixed   statistic    pvalue
## 1 1 157.3255 intensity ~ rough_diagnosis 0.837463338 0.3601233
## 2 2 171.3083 intensity ~ rough_diagnosis 0.090390044 0.7636819
## 3 3 179.1700 intensity ~ rough_diagnosis 0.087791357 0.7670040
## 4 4 181.1523 intensity ~ rough_diagnosis 0.005578811 0.9404602
## 5 5 201.2486 intensity ~ rough_diagnosis 1.846483806 0.1741929
## 6 6 215.2422 intensity ~ rough_diagnosis 0.800201667 0.3710331
## ... and 48 more results

2.3.2 Interpreting the results

Again, we can use topFeatures() to find differentially abundant peaks.

rcc_stest_top <- topFeatures(rcc_stest)

subset(rcc_stest_top, fdr < 0.05)
## DataFrame with 2 rows and 5 columns
##           i        mz statistic      pvalue        fdr
##   <integer> <numeric> <numeric>   <numeric>  <numeric>
## 1        51   884.401   14.2910 0.000156612 0.00675664
## 2        52   885.436   13.4103 0.000250246 0.00675664

This time we find 2 differentially abundant peaks (though one is likely an isotope of the other).

plot(rcc_stest, i=c("m/z = 884.40"=51, "m/z = 885.44"=52), col=dpal("Set1"))

Their mean intensities are significantly higher in tumor tissue.

image(rcc_peaks, mz=885.43, layout=c(2,4), free="xy",
    smooth="bilateral", enhance="adaptive", scale=TRUE)

3 Session information

sessionInfo()
## R version 4.4.0 beta (2024-04-15 r86425)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 22.04.4 LTS
## 
## Matrix products: default
## BLAS:   /home/biocbuild/bbs-3.19-bioc/R/lib/libRblas.so 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_GB              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       
## 
## time zone: America/New_York
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
## [1] Cardinal_3.6.0      S4Vectors_0.42.0    BiocParallel_1.38.0
## [4] BiocGenerics_0.50.0 ProtGenerics_1.36.0 BiocStyle_2.32.0   
## 
## loaded via a namespace (and not attached):
##  [1] sass_0.4.9               tiff_0.1-12              bitops_1.0-7            
##  [4] jpeg_0.1-10              lattice_0.22-6           magrittr_2.0.3          
##  [7] digest_0.6.35            evaluate_0.23            grid_4.4.0              
## [10] bookdown_0.39            fftwtools_0.9-11         fastmap_1.1.1           
## [13] jsonlite_1.8.8           Matrix_1.7-0             ontologyIndex_2.12      
## [16] CardinalWorkflows_1.36.0 matter_2.6.0             DBI_1.2.2               
## [19] biglm_0.9-2.1            tinytex_0.50             BiocManager_1.30.22     
## [22] codetools_0.2-20         jquerylib_0.1.4          abind_1.4-5             
## [25] cli_3.6.2                rlang_1.1.3              CardinalIO_1.2.0        
## [28] Biobase_2.64.0           EBImage_4.46.0           cachem_1.0.8            
## [31] yaml_2.3.8               tools_4.4.0              parallel_4.4.0          
## [34] locfit_1.5-9.9           R6_2.5.1                 png_0.1-8               
## [37] lifecycle_1.0.4          magick_2.8.3             htmlwidgets_1.6.4       
## [40] irlba_2.3.5.1            bslib_0.7.0              Rcpp_1.0.12             
## [43] xfun_0.43                highr_0.10               knitr_1.46              
## [46] htmltools_0.5.8.1        nlme_3.1-164             rmarkdown_2.26          
## [49] compiler_4.4.0           RCurl_1.98-1.14