Class comparison for MSI experiments with Cardinal
Revised: September 12, 2024
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: FALSEWe 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_20251024110408
## mass range: 151.2379 to 994.4166
## centroided: TRUERather 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_20251024110408
## mass range: 157.3255 to 887.4361
## centroided: TRUEThis 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)