In kuwisdelu/matter: Out-of-core statistical computing and signal processing

BiocStyle::markdown()

library(matter)
register(SerialParam())

Introduction

Matter 2 provides a variety of signal processing tools for both uniformly-sampled and nonuniformly-sampled signals in both 1D and 2D, as well as a limited selection of processing tools for signals of arbitrary dimension.

While the primary motivation for implementing these signal processing tools is for use with mass spectrometry imaging experiments, most of the provided functions are broadly applicable to any signal processing domain.

Vocabulary

Matter uses a few consistent terms across its available functions when referring to certain aspects of signal processing. These terms will frequently appear as parameters to signal processing functions.

Dimensionality and domain

Signals may have multiple dimensions. Each dimension corresponds to a domain.

Here is a simple 1-dimensional signal that resembles a mass spectrum:

set.seed(1)
s <- simspec(1)

plot_signal(s, xlab="Index", ylab="Intensity")

A 1-dimensional signal has a single domain.

For a mass spectrum, the domain values are the mass-to-charge ratios, (or m/z-values).

plot_signal(domain(s), s, xlab="m/z", ylab="Intensity")

A 2-dimensional signal has two domains.

For example, for an image (which is a 2D signal), the domain values are the x and y locations of the pixels.

plot_image(volcano)

The dimensionality of a signal is the number of domains.

For example, a mass spectrum collected in tandem with liquid chromatography and ion mobility spectrometry (i.e., LC-IMS-MS) is a 3-dimensional signal. The domains for LC-IMS-MS spectra would include (1) retension times, (2) drift times, and (3) m/z-values.

Note that the values of a signal (e.g., intensities) are not considered as a separate dimension or domain.

Index and domain

In some cases, it is necessary for Matter to distinguish between the observed domain values for a signal and the effective domain values for a signal.

This is most obvious when we need to represent multiple nonuniformly-sampled signals (that may have been observed at different domain values) with a single set of domain values.

set.seed(1)
s1 <- simspec(1)
s2 <- simspec(1)

# spectra with different m/z-values
head(domain(s1))
head(domain(s2))

# create a shared vector of m/z-values
mzr <- range(domain(s1), domain(s2))
mz <- seq(from=mzr[1], to=mzr[2], by=0.2)

# create representations with the same m/z-values
s1 <- sparse_vec(s1, index=domain(s1), domain=mz)
s2 <- sparse_vec(s2, index=domain(s2), domain=mz)

In cases like this, Matter refers to the observed domain values as the signal index, and the effective domain values as the signal domain.

Filtering and smoothing

Smoothing in 1D

Matter provides a family of 1D smoothing methods that follow the naming scheme filt1_*:

• filt1_ma performs moving average filtering

• filt1_conv performs convolutional filtering with custom weights

• filt1_gauss performs Gaussian filtering

• filt1_bi performs bilateral filtering

• filt1_adapt performs adaptive bilateral filtering

• filt1_diff performs nonlinear diffusion filtering

• filt1_guide performs guided filtering

• filt1_pag performs peak-aware guided filtering

• filt1_sg performs Savitzky-Golay filtering

We demonstrate the smoothing results on a simulate signal below:

set.seed(1)
s <- simspec(1, sdnoise=0.25, resolution=500)

p1 <- plot_signal(s)
p2 <- plot_signal(filt1_ma(s))
p3 <- plot_signal(filt1_gauss(s))
p4 <- plot_signal(filt1_bi(s))
p5 <- plot_signal(filt1_diff(s))
p6 <- plot_signal(filt1_guide(s))
p7 <- plot_signal(filt1_pag(s))
p8 <- plot_signal(filt1_sg(s))

plt <- as_facets(list(
    "Original"=p1,
    "Moving average"=p2,
    "Gaussian"=p3,
    "Bilateral"=p4,
    "Diffusion"=p5,
    "Guided"=p6,
    "Peak-aware"=p7,
    "Savitsky-Golay"=p8), nrow=4, ncol=2)

plot(plt)

Now we zoom in on the x-axis to better see the differences in the smoothing.

plot(plt, xlim=c(5800, 6100))

Smoothing in 2D

Matter also provides a family of 2D image smoothing methods that follow the naming scheme filt2_*:

• filt2_ma performs moving average filtering

• filt2_conv performs convolutional filtering with custom weights

• filt2_gauss performs Gaussian filtering

• filt2_bi performs bilateral filtering

• filt2_adapt performs adaptive bilateral filtering

• filt2_diff performs nonlinear diffusion filtering

• filt2_guide performs guided filtering

If a multichannel image is provided (e.g., a 3D array), then these functions will smooth each channel independently.

We demonstrate the smoothing results on the volcano dataset with added noise:

set.seed(1)
img <- volcano + rnorm(length(volcano), sd=2.5)

p1 <- plot_image(img)
p2 <- plot_image(filt2_ma(img))
p3 <- plot_image(filt2_gauss(img))
p4 <- plot_image(filt2_bi(img))
p5 <- plot_image(filt2_diff(img))
p6 <- plot_image(filt2_guide(img))

plt <- as_facets(list(
    "Original"=p1,
    "Moving average"=p2,
    "Gaussian"=p3,
    "Bilateral"=p4,
    "Diffusion"=p5,
    "Guided"=p6), nrow=3, ncol=2)

plot(plt)

One of the major differences between the various smoothing methods is how well they preserve sharp edges.

We use some simple simulated data to demonstrate this:

set.seed(1)
dm <- c(64, 64)
img <- array(rnorm(prod(dm)), dim=dm)
i <- (dm[1] %/% 3):(2 * dm[1] %/% 3)
j <- (dm[2] %/% 3):(2 * dm[2] %/% 3)
img[i,] <- img[i,] + 2
img[,j] <- img[,j] + 2

p1 <- plot_image(img)
p2 <- plot_image(filt2_ma(img))
p3 <- plot_image(filt2_gauss(img))
p4 <- plot_image(filt2_bi(img))
p5 <- plot_image(filt2_diff(img))
p6 <- plot_image(filt2_guide(img))

plt <- as_facets(list(
    "Original"=p1,
    "Moving average"=p2,
    "Gaussian"=p3,
    "Bilateral"=p4,
    "Diffusion"=p5,
    "Guided"=p6), nrow=3, ncol=2)

plot(plt)

Smoothing using KNN

A limited set of smoothing filters are also available for N-dimensional signals that have been nonuniformly sampled. For example, this could include spatial signals at arbitrary (non-gridded) locations or mass spectra with an ion mobility dimension.

These follow the naming scheme filtn_* and include:

• filtn_ma performs moving average filtering

• filtn_conv performs convolutional filtering with custom weights

• filtn_gauss performs Gaussian filtering

• filtn_bi performs bilateral filtering

• filtn_adapt performs adaptive bilateral filtering

Because these smoothing functions rely on k-nearest-neighbor search to determine the smoothing windows, they are more computationally intensive than the standard 1D and 2D filters.

Contrast enhancement

Contrast enhancement improves the contrast of an image. This can correct for hotspots caused by multiplicative variance.

To demonstrate this, we will add multiplicative noise to the volcano dataset:

set.seed(1)
img <- volcano + rlnorm(length(volcano), sd=1.5)

Contrast enhancement functions follow the naming scheme enhance_* and includes:

• enhance_adj adjusts contrast by clamping extreme values

• enhance_hist performs histogram equalization

• enhance_adapt performs contrast-limited adaptive histogram equalization (CLAHE)

p1 <- plot_image(img)
p2 <- plot_image(enhance_adj(img))
p3 <- plot_image(enhance_hist(img))
p4 <- plot_image(enhance_adapt(img))

plt <- as_facets(list(
    "Original"=p1,
    "Adjust"=p2,
    "Histogram"=p3,
    "CLAHE"=p4), nrow=2, ncol=2)

plot(plt, scale=TRUE)

Using contrast enhancement can dramatically improve the visual interpretation of the image.

Note that the contrast-limited adaptive histogram equalization (CLAHE) method is most useful when the required contrast correction is very different across different regions of the image, which is not the case here.

Rescaling (normalization)

Rescaling the signal is often necessary for various normalization methods.

Rescaling functions follow the naming scheme rescale_* and include:

• rescale_rms rescales based on the root-mean-square of the signal

• rescale_sum rescales based on the sum of (the absolute values of) the signal

• rescale_ref rescales according to a specific sample (i.e., a reference value at a specific index)

• rescale_range rescales the lower and upper limits of the signal values

• rescale_iqr rescales the lower and upper limits of the signal values

As they only rescale the signal we won't visualize these methods here.

Continuum estimation

Continuum estimation is most commonly used to estimate a signal's baseline.

set.seed(1)
s <- simspec(1, baseline=5, resolution=500)

Baseline estimation functions follow the naming scheme estbase_* and include:

• estbase_loc interpolates the continuum from local extrema

• estbase_hull estimates an upper or lower convex hull

• estbase_snip performs sensitive nonlinear iterative peak clipping (SNIP)

• estbase_med estimates the continuum form running medians

We demonstrate baseline removal below:

p1 <- plot_signal(s)
p2 <- plot_signal(s - estbase_loc(s))
p3 <- plot_signal(s - estbase_hull(s))
p4 <- plot_signal(s - estbase_snip(s))

plt <- as_facets(list(
    "Original"=p1,
    "Local minima"=p2,
    "Convex hull"=p3,
    "SNIP"=p4), nrow=2, ncol=2)

plot(plt)

Warping and alignment

It is often necessary to warp signals to align their features (e.g., peaks) in the presence of variance in their observed domain values.

This is necessary, for example, to recalibrate mass spectra with a large amount of mass error, or to co-register images across different image modalities.

Warping in 1D

Matter provides a few signal warping functions that follow the naming scheme warp1_*:

• warp1_loc performs warping based on local extrema

• warp1_dtw performs dynamic time warping

• warp1_cow performs correlation optimized warping

These functions can be align two signals as shown below.

First, we need to simulate some signals in need of alignment:

set.seed(1)
s <- simspec(8, sdx=5e-4)

plot_signal(domain(s), s, by=NULL, group=1:8,
    xlim=c(1250, 1450))

Next, we need to generate a reference signal to use when aligning all of the signals. We do this by calculating the mean signal and using it as the reference:

ref <- rowMeans(s)
s2 <- apply(s, 2L, warp1_loc,
    y=ref, events="max", tol=2e-3, tol.ref="y")

plot_signal(domain(s), s2, by=NULL, group=1:8,
    xlim=c(1250, 1450))

Note that signal warping can be quite sensitive to the tolerance (i.e., how much the signal is allowed to shift in either direction).

Warping in 2D

A warping function for 2D images is also available.

First, we need to generate images in need of co-registration:

img1 <- trans2d(volcano, rotate=15, translate=c(-5, 5))

plot_image(list(volcano, img1))

Now, we can warp the transformed image to align it with the original image:

img2 <- warp2_trans(img1, volcano)

plot_image(list(volcano, img2))

Note that this is method is currently limited to affine transformations.

Peak processing

Matter provides a family of functions for detecting peaks in a signal and processing peak lists from multiple signals.

To demonstrate peak processing, we begin by simulating a signal.

set.seed(1)
s <- simspec(1)

Local maxima

A typical first step in peak processing is detection of local maxima to select peak candidates.

Local maxima in 1D

The locmax and locmin functions perform detection of local extreme based on a sliding window approach.

A sample is considered a local maximum if it is greater than all elements to its left within the window and greater than or equal to all elements in the window to its right within the window.

The default window has width=5:

i <- locmax(s)

plot_signal(domain(s), s, xlim=c(900, 1100))
lines(domain(s)[i], s[i], type="h", col="red")

Below we use a wider window with width=15.

i <- locmax(s, width=15)

plot_signal(domain(s), s, xlim=c(900, 1100))
lines(domain(s)[i], s[i], type="h", col="red")

In a noisy signal, there will be many more local maxima than true peaks.

Local maxima using KNN

For signals with multiple dimensions, the knnmax and knnmin functions perform detection of local extrema using k-nearest-neighbors.

Below, we detect the local maxima in a simulated dataset:

set.seed(1)
dm <- c(64, 64)
img <- array(rnorm(prod(dm)), dim=dm)
w <- 100 * dnorm(seq(-3, 3)) %o% dnorm(seq(-3, 3))
img[1:7,1:7] <- img[1:7,1:7] + w
img[11:17,11:17] <- img[11:17,11:17] + w
img[21:27,21:27] <- img[21:27,21:27] + w
img[21:27,41:47] <- img[21:27,41:47] + w
img[51:57,31:37] <- img[51:57,31:37] + w
img[41:47,51:57] <- img[41:47,51:57] + w

plot_image(img)

coord <- expand.grid(lapply(dim(img), seq_len))

i <- knnmax(img, index=coord, k=49)
i <- which(i, arr.ind=TRUE)

plot_image(img)
points(i[,1], i[,2], pch=4, lwd=2, cex=2, col="red")

As in 1-dimension, there are many more local maxima than true peaks, but the true peaks (i.e., the bright spots) are detected.

Noise estimation

To distinguish true peaks from false peaks, we need to estimate the noise level in the signal.

Matter provides a variety of noise estimation functions that follow the naming scheme estnoise_*:

• estnoise_diff estimates noise from the differences between the signal and its derivative

• estnoise_filt uses dynamic filtering-based to distinguish between noise peaks and true peaks

• estnoise_sd estimates noise from the standard deviation of the difference between a smoothed signal and the original signal

• estnoise_mad estimates noise from the mean absolute deviation of the difference between a smoothed signal and the original signal

• estnoise_quant estimates noise from a quantile of the difference between a smoothed signal and the original signal

p1 <- plot_signal(domain(s), s,
    xlim=c(900, 1100), ylim=c(0, 15))

p2 <- add_mark(p1, "lines",
    x=domain(s), y=estnoise_diff(s),
    params=list(color="blue", linewidth=2))

p3 <- add_mark(p1, "lines",
    x=domain(s), y=estnoise_filt(s),
    params=list(color="blue", linewidth=2))

p4 <- add_mark(p1, "lines",
    x=domain(s), y=estnoise_quant(s),
    params=list(color="blue", linewidth=2))

p5 <- add_mark(p1, "lines",
    x=domain(s), y=estnoise_sd(s),
    params=list(color="blue", linewidth=2))

p6 <- add_mark(p1, "lines",
    x=domain(s), y=estnoise_mad(s),
    params=list(color="blue", linewidth=2))

as_facets(list(
    "Original"=p1,
    "Difference"=p2,
    "Filtering"=p3,
    "Quantile"=p4,
    "SD"=p5,
    "MAD"=p6), nrow=3, ncol=2)

By default, a constant noise level is estimated for the entire signal domain. A variable noise level can be estimated for different regions of the signal by specifying the nbins argument.

p1 <- plot_signal(domain(s), s,
    xlim=c(900, 1100), ylim=c(0, 15))

p2 <- add_mark(p1, "lines",
    x=domain(s), y=estnoise_diff(s, nbins=20),
    params=list(color="blue", linewidth=2))

p3 <- add_mark(p1, "lines",
    x=domain(s), y=estnoise_filt(s, nbins=20),
    params=list(color="blue", linewidth=2))

p4 <- add_mark(p1, "lines",
    x=domain(s), y=estnoise_quant(s, nbins=20),
    params=list(color="blue", linewidth=2))

p5 <- add_mark(p1, "lines",
    x=domain(s), y=estnoise_sd(s, nbins=20),
    params=list(color="blue", linewidth=2))

p6 <- add_mark(p1, "lines",
    x=domain(s), y=estnoise_mad(s, nbins=20),
    params=list(color="blue", linewidth=2))

as_facets(list(
    "Original"=p1,
    "Difference"=p2,
    "Filtering"=p3,
    "Quantile"=p4,
    "SD"=p5,
    "MAD"=p6), nrow=3, ncol=2)

Peak detection

The findpeaks function streamlines the process of peak detection. It detects local maxima using locmax, optionally estimates noise using a specified estnoise_* function, and optionally filters the peaks based on a signal-to-noise ratio threshold set via snr.

i1 <- findpeaks(s, snr=3, noise="diff")

p1a <- plot_signal(domain(s), s, group="Signal")
p1b <- plot_signal(domain(s)[i1], s[i1], group="Peaks",
    isPeaks=TRUE, annPeaks="circle")

p1 <- as_layers(p1a, p1b)

i2 <- findpeaks(s, snr=3, noise="filt")

p2a <- plot_signal(domain(s), s, group="Signal")
p2b <- plot_signal(domain(s)[i2], s[i2], group="Peaks",
    isPeaks=TRUE, annPeaks="circle")

p2 <- as_layers(p2a, p2b)

i3 <- findpeaks(s, snr=3, noise="sd")

p3a <- plot_signal(domain(s), s, group="Signal")
p3b <- plot_signal(domain(s)[i3], s[i3], group="Peaks",
    isPeaks=TRUE, annPeaks="circle")

p3 <- as_layers(p3a, p3b)

i4 <- findpeaks(s, snr=3, noise="mad")

p4a <- plot_signal(domain(s), s, group="Signal")
p4b <- plot_signal(domain(s)[i4], s[i4], group="Peaks",
    isPeaks=TRUE, annPeaks="circle")

p4 <- as_layers(p4a, p4b)

plt <- as_facets(list(
    "Difference"=p1,
    "Filtering"=p2,
    "SD"=p3,
    "MAD"=p4))

plot(plt, xlim=c(1200, 1600))

Different noise estimation methods may have different sensitivity and specificity when detecting peaks in different kinds of signals.

Peak binning

After detecting peaks from multiple signals, it is often necessary to determine whether two (or more) peaks correspond to the same feature or not. This is typically done by binning the peaks.

For example, in the context of mass spectrometry, there is typically some small amount of error in the observed m/z-values of the peaks. Peak binning is necessary to know which peaks from different spectra refer to the same molecules.

set.seed(1)
s <- simspec(4)

peaklist <- apply(s, 2, findpeaks, snr=3)

peaks <- binpeaks(peaklist)
as.vector(peaks)

The peaks are binned to a domain that is automatically generated based on the range of the detected peaks (if domain is not explicitly specified).

The bins of the domain are spaced according to a tolerance tol. If tol is not specified, then it is estimated as half of the median of the minimum gap between same-signal peaks.

The value returned by binpeaks above are the binned peaks, as determined by the domain bins where at least one peak was observed. The per-signal peak locations are averaged within each bin to determine the final peak locations.

Peak merging

After peak binning, there may still be some peaks that need to be merged.

Whether this is necessary depends on the resolution of the bins that were used when binning the peaks. Higher resolution bins (preferred for the sake of accuracy) may require merging bins that are likely the same peak.

peaks <- mergepeaks(peaks, tol=1)
as.vector(peaks)

As can be seen above, the peaks that were binned at 1000 and 1001 previously have been merged into a single peak at 1000.75.

Note that the number of observed peaks at each bin are preserved and used to calculate the new peak location, which is why the resulting peak is located at 1000.75 instead of 1000.5.

Session information

sessionInfo()

kuwisdelu/matter documentation built on Oct. 19, 2024, 10:31 a.m.