R/process-peakPick.R

In Cardinal: A mass spectrometry imaging toolbox for statistical analysis

#### Peak picking methods ####
## ---------------------------

setMethod("peakPick", "MSImageSet",
	function(object, method = c("simple", "adaptive", "limpic"),
		...,
		pixel = pixels(object),
		plot = FALSE)
	{
		.Deprecated_Cardinal1()
		if ( centroided(object) )
			.stop("peakPick: Data already centroided. Peak picking will not be performed.")
		fun <- peakPick.method(method)
		.message("peakPeak: Using method = ", match.method(method))
		.time.start()
		peaks <- pixelApply(object, function(s, ...) {
			peakPick.do(s, object, .Index, fun, plot, ...)
		}, .pixel=pixel, ..., .use.names=FALSE, .simplify=FALSE)
		peakData <- lapply(seq_along(pixel), function(i) {
			intensityArray <- iData(object)[peaks[[i]], pixel[i]]
			names(intensityArray) <- featureNames(object)[peaks[[i]]]
			intensityArray
		})
		mzData <- lapply(seq_along(pixel), function(i) {
			mzArray <- mz(object)[peaks[[i]]]
			names(mzArray) <- featureNames(object)[peaks[[i]]]
			mzArray
		})
		object <- object[,pixel]
		peakData <- new("Hashmat", data=peakData, keys=featureNames(object),
			dim=c(length(features(object)), length(peakData)))
		mzData <- new("Hashmat", data=mzData, keys=featureNames(object),
			dim=c(length(features(object)), length(mzData)))
		peakData(imageData(object)) <- peakData
		mzData(imageData(object)) <- mzData
		peakPicking(processingData(object)) <- match.method(method)
		.message("peakPick: Done.")
		.time.stop()
		object
	})

peakPick.do <- function(s, object, pixel, f, plot, ...) {
	pout <- f(s, ...)
	if ( plot ) {
		wrap(plot(object, s ~ mz, pixel=pixel, col="gray",
			ylab="Intensity", strip=FALSE, ...),
			..., signature=f)
		wrap(lines(mz(object), pout$noise, col="blue", ...),
			..., signature=f)
		wrap(lines(mz(object)[pout$peaks], s[pout$peaks], col="red", type='h', ...),
			..., signature=f)
	}
	pout$peaks
}

peakPick.method <- function(method, name.only=FALSE) {
	if ( is.character(method) || is.null(method) ) {
		options <- c("simple", "adaptive", "limpic")
		method <- match.method(method, options)
		if ( name.only )
			return(method)
		method <- switch(method,
			simple = peakPick.simple,
			adaptive = peakPick.adaptive,
			limpic = peakPick.limpic,
			match.fun(method))
	}
	match.fun(method)
}

peakPick.simple <- function(x, SNR=6, window=5, blocks=100, ...) {
	noise <- .estimateNoiseSimple(x, blocks=blocks)
	maxs <- locmax(x, halfWindow=window %/% 2)
	peaks <- intersect(maxs, which(x / noise >= SNR))
	list(peaks=peaks, noise=noise)
}

peakPick.adaptive <- function(x, SNR=6, window=5, blocks=100, spar=1, ...) {
	noise <- .estimateNoiseAdaptive(x, blocks=blocks, spar=spar)
	maxs <- locmax(x, halfWindow=window %/% 2)
	peaks <- intersect(maxs, which(x / noise >= SNR))
	list(peaks=peaks, noise=noise)
}

peakPick.limpic <- function(x, SNR=6, window=5, blocks=100, thresh=0.75, ...) {
	t <- seq_along(x)
	xint <- split_blocks(x, blocks=blocks)
	# identify flat regions of spectrum
	kurt <- sapply(xint, kurtosis) - 3
	kurt[is.nan(kurt)] <- -Inf
	means <- sapply(xint, mean)
	is.flat <- kurt < 1 & means < mean(x)
	# estimate noise
	noiseval <- sapply(xint, sd)[is.flat]
	noiseval <- c(median(noiseval, na.rm=TRUE), noiseval, median(noiseval, na.rm=TRUE))
	noiseidx <- floor(seq(from=blocks/2, to=length(x) - blocks/2,
		length.out=length(xint)))[is.flat]
	noiseidx <- c(1, noiseidx, length(x))
	noise <- interp1(x=t[noiseidx], y=noiseval, xi=t, method="linear")
	# find local maxima
	halfWindow <- floor(window / 2)
	maxs <- locmax(x, halfWindow=window %/% 2)
	peaks <- intersect(maxs, which(x / noise >= SNR))
	# eliminate smooth peaks
	cutoff <- ceiling(halfWindow / 2) * quantile(abs(diff(x)), thresh) / 2
	too.flat <- sapply(peaks, function(i) {
		if ( (i - halfWindow) < 0 || (i + halfWindow) > length(x) ) {
			TRUE	
		} else {
			ratio <- x[i] - 0.5 * (x[i - halfWindow] + x[i + halfWindow])
			ratio < cutoff
		}
	})
	# return peak list
	peaks <- setdiff(peaks, too.flat)
	list(peaks=peaks, noise=noise)
}

# Estimate noise in a signal

.estimateNoiseSimple <- function(x, blocks=100) {
	t <- seq_along(x)
	xint <- split_blocks(x, blocks=blocks)
	kurt <- sapply(xint, kurtosis) - 3
	if ( all(kurt >= 1, na.rm=TRUE) ) {
		noise <- min(sapply(xint, sd), na.rm=TRUE)
		noise <- rep(noise, length(x))
		if ( anyNA(kurt) )
			.warning("kurtosis could not be estimated; ",
				"try method = 'mad' instead")
	} else {
		kurt[is.na(kurt)] <- Inf
		noise <- mean(sapply(xint, sd)[kurt < 1], na.rm=TRUE)
		noise <- rep(noise, length(x))
	}
	noise
}

.estimateNoiseAdaptive <- function(x, blocks=100, spar=1) {
	t <- seq_along(x)
	tint <- split_blocks(t, blocks=blocks)
	xint <- split_blocks(x, blocks=blocks)
	kurt <- sapply(xint, kurtosis) - 3
	if ( all(kurt >= 1, na.rm=TRUE) ) {
		noise <- min(sapply(xint, sd), na.rm=TRUE)
		noise <- rep(noise, length(x))
		if ( anyNA(kurt) )
			.warning("kurtosis could not be estimated; ",
				"try method = 'mad' instead")
	} else {
		kurt[is.na(kurt)] <- Inf
		noiseval <- sapply(xint, sd)[kurt < 1]
		noiseidx <- sapply(tint, mean)[kurt < 1]
		noise <- interp1(noiseidx, noiseval, xi=t, method="linear",
			extrap=mean(noiseval, na.rm=TRUE))
		noise <- smooth.spline(x=t, y=noise, spar=1)$y
	}
	noise
}