| Title: | Peak Detection for Mass Spectrometry data using wavelet-based algorithms |
|---|---|
| Description: | Peak Detection in Mass Spectrometry data is one of the important preprocessing steps. The performance of peak detection affects subsequent processes, including protein identification, profile alignment and biomarker identification. Using Continuous Wavelet Transform (CWT), this package provides a reliable algorithm for peak detection that does not require any type of smoothing or previous baseline correction method, providing more consistent results for different spectra. See <doi:10.1093/bioinformatics/btl355} for further details. |
| Authors: | Pan Du [aut],
Warren Kibbe [aut],
Simon Lin [aut],
Sergio Oller Moreno [aut, cre]
|
| Maintainer: | Sergio Oller Moreno <[email protected]> |
| License: | LGPL (>= 2) |
| Version: | 1.73.1 |
| Built: | 2024-12-18 08:27:46 UTC |
| Source: | https://github.com/zeehio/MassSpecWavelet |
MassSpecWavelet R package is aimed to detect peaks on Mass Spectrometry (MS) data using Continuous Wavelet Transform (CWT).
Pan Du, Simon Lin
Du, P., Kibbe, W.A. and Lin, S.M. (2006) Improved peak detection in mass spectrum by incorporating continuous wavelet transform-based pattern matching, Bioinformatics, 22, 2059-2065.
Useful links:
Report bugs at http://github.com/zeehio/MassSpecWavelet/issues
data(exampleMS) SNR.Th <- 3 peakInfo <- peakDetectionCWT(exampleMS, SNR.Th = SNR.Th) majorPeakInfo <- peakInfo$majorPeakInfo peakIndex <- majorPeakInfo$peakIndex plotPeak(exampleMS, peakIndex, main = paste("Identified peaks with SNR >", SNR.Th))data(exampleMS) SNR.Th <- 3 peakInfo <- peakDetectionCWT(exampleMS, SNR.Th = SNR.Th) majorPeakInfo <- peakInfo$majorPeakInfo peakIndex <- majorPeakInfo$peakIndex plotPeak(exampleMS, peakIndex, main = paste("Identified peaks with SNR >", SNR.Th))
CWT(Continuous Wavelet Transform) with Mexican Hat wavelet (by default) to match the peaks in Mass Spectrometry spectrum
cwt(ms, scales = 1, wavelet = "mexh")cwt(ms, scales = 1, wavelet = "mexh")
ms |
Mass Spectrometry spectrum (a vector of MS intensities) |
scales |
a vector represents the scales at which to perform CWT. See
the Details section. Additionally, a |
wavelet |
The wavelet base, Mexican Hat by default. User can provide
wavelet |
The default mother wavelet is a Mexican Hat evaluated in the range
using 1024 points (a step of roughly 1/64):
The of the mother Mexican Hat is of one x unit.
The scaled wavelet is a downsampled version of the mother wavelet. The
scale determines how many samples per unit are taken. For
instance, with the default Mexican Hat wavelet, a scales = 2 will
evaluate the range sampling twice per unit, this is with a
step of 0.5. This generates a 33 points long scaled wavelet. Choosing this
type of scaling is convenient because the scaled wavelet becomes a wavelet
of points. Using the default wavelet, a
scales smaller than 1 will show sampling issues, while a
scales larger than 64 will resample points from the original mother
wavelet. If you need to use scales larger than 64, consider providing
your own mother wavelet. See the examples.
According to doi:10.1063/1.3505103, if
your spectrum has a gaussian peak shape of variance points then
the scales range should cover . If your spectrum has a
Lorentzian peak shape of half-width-half-maximum points then the
scales range should cover . They also suggest using a
spacing. Take these values as suggestions for your
analysis.
The return is the 2-D CWT coefficient matrix, with column names as the scale. Each column is the CWT coefficients at that scale. If the scales are too big for the given signal, the returned matrix may include less columns than the given scales.
Pan Du, Simon Lin
data(exampleMS) scales <- seq(1, 64, 3) wCoefs <- cwt(exampleMS[5000:11000], scales = scales, wavelet = "mexh") ## Plot the 2-D CWT coefficients as image (It may take a while!) xTickInterval <- 1000 image(5000:11000, scales, wCoefs, col = terrain.colors(256), axes = FALSE, xlab = "m/z index", ylab = "CWT coefficient scale", main = "CWT coefficients" ) axis(1, at = seq(5000, 11000, by = xTickInterval)) axis(2, at = c(1, seq(10, 64, by = 10))) box() ## Provide a larger wavelet: scales <- c(seq(1, 64, 3), seq(72, 128, 8)) psi_xval <- seq(-8, 8, length.out = 4096) psi <- (2 / sqrt(3) * pi^(-0.25)) * (1 - psi_xval^2) * exp(-psi_xval^2 / 2) wCoefs <- cwt(exampleMS[5000:11000], scales = scales, wavelet = rbind(psi_xval, psi)) xTickInterval <- 1000 image(5000:11000, scales, wCoefs, col = terrain.colors(256), axes = FALSE, xlab = "m/z index", ylab = "CWT coefficient scale", main = "CWT coefficients" ) axis(1, at = seq(5000, 11000, by = xTickInterval)) axis(2, at = c(1, seq(10, 128, by = 10))) box() ## Custom log1.18 spaced scales: get_scales <- function(scale_min, scale_max, num_scales) { (seq(0, 1, length.out = num_scales)^1.18) * (scale_max - scale_min) + scale_min } scales <- get_scales(scale_min = 1, scale_max = 64, num_scales = 32) print(scales) # For detecting a gaussian peak of 10 points: gaussian_peak_sigma <- 10 # points scales <- get_scales(scale_min = 1, scale_max = 1.9 * gaussian_peak_sigma, num_scales = 32) print(scales) # For detecting a lorentzian peak of 10 points: lorentzian_peak_gamma <- 10 # points scales <- get_scales(scale_min = 1, scale_max = 2.9 * lorentzian_peak_gamma, num_scales = 32) print(scales)data(exampleMS) scales <- seq(1, 64, 3) wCoefs <- cwt(exampleMS[5000:11000], scales = scales, wavelet = "mexh") ## Plot the 2-D CWT coefficients as image (It may take a while!) xTickInterval <- 1000 image(5000:11000, scales, wCoefs, col = terrain.colors(256), axes = FALSE, xlab = "m/z index", ylab = "CWT coefficient scale", main = "CWT coefficients" ) axis(1, at = seq(5000, 11000, by = xTickInterval)) axis(2, at = c(1, seq(10, 64, by = 10))) box() ## Provide a larger wavelet: scales <- c(seq(1, 64, 3), seq(72, 128, 8)) psi_xval <- seq(-8, 8, length.out = 4096) psi <- (2 / sqrt(3) * pi^(-0.25)) * (1 - psi_xval^2) * exp(-psi_xval^2 / 2) wCoefs <- cwt(exampleMS[5000:11000], scales = scales, wavelet = rbind(psi_xval, psi)) xTickInterval <- 1000 image(5000:11000, scales, wCoefs, col = terrain.colors(256), axes = FALSE, xlab = "m/z index", ylab = "CWT coefficient scale", main = "CWT coefficients" ) axis(1, at = seq(5000, 11000, by = xTickInterval)) axis(2, at = c(1, seq(10, 128, by = 10))) box() ## Custom log1.18 spaced scales: get_scales <- function(scale_min, scale_max, num_scales) { (seq(0, 1, length.out = num_scales)^1.18) * (scale_max - scale_min) + scale_min } scales <- get_scales(scale_min = 1, scale_max = 64, num_scales = 32) print(scales) # For detecting a gaussian peak of 10 points: gaussian_peak_sigma <- 10 # points scales <- get_scales(scale_min = 1, scale_max = 1.9 * gaussian_peak_sigma, num_scales = 32) print(scales) # For detecting a lorentzian peak of 10 points: lorentzian_peak_gamma <- 10 # points scales <- get_scales(scale_min = 1, scale_max = 2.9 * lorentzian_peak_gamma, num_scales = 32) print(scales)
An example mass spectrum from CAMDA 2006. All-in-1 Protein Standard II (Ciphergen Cat. \# C100-0007) were measured on Ciphergen NP20 chips. There are 7 polypeptides in the sample with m/z values of 7034, 12230, 16951, 29023, 46671, 66433, 147300.
A numeric vector represents the mass spectrum with equal sample intervals.
CAMDA, CAMDA 2006 Competition Data Set. 2006, http://camda.duke.edu.
Compared to the rest of the package, this is a rather experimental function. If you plan to use it or are interested in it, please open an issue at https://github.com/zeehio/MassSpecWavelet/issues to show your interest.
findLocalMaxWinSize(x, capWinSize = NA)findLocalMaxWinSize(x, capWinSize = NA)
x |
A numeric vector. |
capWinSize |
the maximum window size to report. |
An integer vector y of the same length as x. y[i] will be the
size of the largest window on x containing x[i] where:
x[i] is a local maximum or a center of a plateau
x[i] is not at a window border
Optionally, if capWinSize is a positive integer, the maximum window size
is capped to that value, to increase performance. Use this in case you just
want to check if there exists a window of that size.
@export
@examples
x <- c(1, 2, 3, 2, 1)
findLocalMaxWinSize(x)
Identify the local maximum of each column in 2-D CWT coefficients matrix by using a slide window. The size of slide window linearly changes from the coarse scale (bigger window size) to detail scale. The scale of CWT increases with the column index.
getLocalMaximumCWT( wCoefs, minWinSize = 5, amp.Th = 0, isAmpThreshRelative = FALSE, exclude0scaleAmpThresh = FALSE )getLocalMaximumCWT( wCoefs, minWinSize = 5, amp.Th = 0, isAmpThreshRelative = FALSE, exclude0scaleAmpThresh = FALSE )
wCoefs |
2-D CWT coefficients, each column corresponding to CWT coefficient at one scale. The column name is the scale. |
minWinSize |
The minimum slide window size used. |
amp.Th |
The minimum peak amplitude. |
isAmpThreshRelative |
Whether |
exclude0scaleAmpThresh |
When computing the relative |
return a matrix with same dimension as CWT coefficient matrix, wCoefs. The local maxima are marked as 1, others are 0.
Pan Du
data(exampleMS) scales <- seq(1, 64, 3) wCoefs <- cwt(exampleMS[5000:11000], scales = scales, wavelet = "mexh") localMax <- getLocalMaximumCWT(wCoefs) plotLocalMax(localMax)data(exampleMS) scales <- seq(1, 64, 3) wCoefs <- cwt(exampleMS[5000:11000], scales = scales, wavelet = "mexh") localMax <- getLocalMaximumCWT(wCoefs) plotLocalMax(localMax)
Identify ridges by connecting the local maximum of 2-D CWT coefficients from
the coarse scale to detail scale. The local maximum matrix is returned from
getLocalMaximumCWT()
getRidge( localMax, iInit = ncol(localMax), step = -1, iFinal = 1, minWinSize = 5, gapTh = 3, skip = NULL, scaleToWinSize = "doubleodd" )getRidge( localMax, iInit = ncol(localMax), step = -1, iFinal = 1, minWinSize = 5, gapTh = 3, skip = NULL, scaleToWinSize = "doubleodd" )
localMax |
The local maximum matrix is returned from
|
iInit |
The start column to search ridge. By default, it starts from the coarsest scale level. |
step |
Search step. -1 by default, which means searching from coarse scale to detail scale column by column. |
iFinal |
The final column index of search ridge. |
minWinSize |
The minimum slide window size used. |
gapTh |
The gap allowed during searching for ridge. 3 by default. |
skip |
The column to be skipped during search. |
scaleToWinSize |
How scales should be mapped to window sizes. Traditionally,
MassSpecWavelet used the |
Return a list of ridge. As some ridges may end at the scale larger than 1, in order to keep the uniqueness of the ridge names, we combined the smallest scale of the ridge and m/z index of the peak at that scale together to name the ridges. For example the ridge name "1\_653" means the peak ridge ends at the CWT scale 1 with m/z index 653 at scale 1.
Pan Du, Simon Lin
Du, P., Kibbe, W.A. and Lin, S.M. (2006) Improved peak detection in mass spectrum by incorporating continuous wavelet transform-based pattern matching, Bioinformatics, 22, 2059-2065.
getLocalMaximumCWT(), identifyMajorPeaks()
data(exampleMS) scales <- seq(1, 64, 3) wCoefs <- cwt(exampleMS[5000:11000], scales = scales, wavelet = "mexh") localMax <- getLocalMaximumCWT(wCoefs) ridgeList <- getRidge(localMax) plotRidgeList(ridgeList)data(exampleMS) scales <- seq(1, 64, 3) wCoefs <- cwt(exampleMS[5000:11000], scales = scales, wavelet = "mexh") localMax <- getLocalMaximumCWT(wCoefs) ridgeList <- getRidge(localMax) plotRidgeList(ridgeList)
Estimate the length of the ridge line, which is composed of local maxima at adjacent CWT scales. The ridge line is cut off at the end point, whose amplitude divided by the maximum ridge amplitude is larger than the cutoff amplitude ratio threshold (0.5 by default).
getRidgeLength(ridgeList, Th = 0.5)getRidgeLength(ridgeList, Th = 0.5)
ridgeList |
a list of identified ridges |
Th |
the cutoff amplitude ratio (the amplitude divided by the maximum amplitude of the ridge) threshold of the ridge line end. |
a vector of estimated ridge length
Pan Du
stopifnot(getRidgeLength(list(c(5,4,3,2,1), c(5,3,1))) == c(3,2))stopifnot(getRidgeLength(list(c(5,4,3,2,1), c(5,3,1))) == c(3,2))
Get the CWT coefficient values corresponding to the peak ridge
getRidgeValue(ridgeList, wCoefs, skip = 0)getRidgeValue(ridgeList, wCoefs, skip = 0)
ridgeList |
a list of ridge lines |
wCoefs |
2-D CWT coefficients |
skip |
the CWT scale level to be skipped, by default the 0 scale level (raw spectrum) is skipped. |
A list of ridge values corresponding to the input ridgeList.
Pan Du
Indentify the peaks based on the ridge list (returned by
getRidge()) in 2-D CWT coefficient matrix and estimated Signal
to Noise Ratio (SNR)
identifyMajorPeaks( ms, ridgeList, wCoefs, scales = as.numeric(colnames(wCoefs)), SNR.Th = 3, peakScaleRange = 5, ridgeLength = 32, nearbyPeak = FALSE, nearbyWinSize = ifelse(nearbyPeak, 150, 100), winSize.noise = 500, SNR.method = "quantile", minNoiseLevel = 0.001, excludeBoundariesSize = nearbyWinSize/2 )identifyMajorPeaks( ms, ridgeList, wCoefs, scales = as.numeric(colnames(wCoefs)), SNR.Th = 3, peakScaleRange = 5, ridgeLength = 32, nearbyPeak = FALSE, nearbyWinSize = ifelse(nearbyPeak, 150, 100), winSize.noise = 500, SNR.method = "quantile", minNoiseLevel = 0.001, excludeBoundariesSize = nearbyWinSize/2 )
ms |
the mass spectrometry spectrum |
ridgeList |
returned by |
wCoefs |
2-D CWT coefficients |
scales |
scales of CWT, by default it is the colnames of wCoefs |
SNR.Th |
threshold of SNR |
peakScaleRange |
the CWT scale range of the peak. |
ridgeLength |
the maximum ridge scale of the major peaks. |
nearbyPeak |
determine whether to include the small peaks close to large major peaks |
nearbyWinSize |
the window size to determine the nearby peaks. Only effective when nearbyPeak is true. |
winSize.noise |
the local window size to estimate the noise level. |
SNR.method |
method to estimate noise level. Currently, only 95 percentage quantile is supported. |
minNoiseLevel |
the minimum noise level used in calculating SNR, i.e., if the estimated noise level is less than "minNoiseLevel", it will use "minNoiseLevel" instead. If the noise level is less than 0.5, it will be treated as the ratio to the maximum amplitude of the spectrum. |
excludeBoundariesSize |
number of points at each boundary of the ms signal that will be excluded in search for peaks to avoid boundary effects. |
The determination of the peaks is based on three rules: Rule 1: The maximum ridge scale of the peak should larger than a certain threshold Rule 1.1: Based on the scale of the peak (corresponding to the maximum value of the peak ridge) should be within certain range Rule 2: Based on the peak SNR Rule 3: The peak should not appear at the boundaries of the signal.
Return a list with following elements:
peakIndex |
the m/z indexes of the identified peaks |
peakCenterIndex |
the m/z indexes of peak centers, which correspond to the maximum on the ridge. peakCenterIndex includes all the peaks, not just the identified major peaks. |
peakCenterValue |
the CWT coefficients (the maximum on the ridge) corresponding to peakCenterIndex |
peakSNR |
the SNR of the peak, which is the ratio of peakCenterValue and noise level |
peakScale |
the estimated scale of the peak, which corresponds to the peakCenerIndex |
potentialPeakIndex |
the m/z indexes of all potential peaks, which satisfy all requirements of a peak without considering its SNR. Useful, if you want to change to a lower SNR threshold later. |
allPeakIndex |
the m/z indexes of all the peaks, whose order is the same as peakCenterIndex, peakCenterValue, peakSNR and peakScale. |
All of these return elements have peak names, which are the same as the
corresponding peak ridges. see getRidge() for details.
Pan Du, Simon Lin
Du, P., Kibbe, W.A. and Lin, S.M. (2006) Improved peak detection in mass spectrum by incorporating continuous wavelet transform-based pattern matching, Bioinformatics, 22, 2059-2065.
peakDetectionCWT(), tuneInPeakInfo()
data(exampleMS) scales <- seq(1, 64, 3) wCoefs <- cwt(exampleMS, scales = scales, wavelet = "mexh") localMax <- getLocalMaximumCWT(wCoefs) ridgeList <- getRidge(localMax) SNR.Th <- 3 majorPeakInfo <- identifyMajorPeaks(exampleMS, ridgeList, wCoefs, SNR.Th = SNR.Th) ## Plot the identified peaks peakIndex <- majorPeakInfo$peakIndex plotPeak(exampleMS, peakIndex, main = paste("Identified peaks with SNR >", SNR.Th))data(exampleMS) scales <- seq(1, 64, 3) wCoefs <- cwt(exampleMS, scales = scales, wavelet = "mexh") localMax <- getLocalMaximumCWT(wCoefs) ridgeList <- getRidge(localMax) SNR.Th <- 3 majorPeakInfo <- identifyMajorPeaks(exampleMS, ridgeList, wCoefs, SNR.Th = SNR.Th) ## Plot the identified peaks peakIndex <- majorPeakInfo$peakIndex plotPeak(exampleMS, peakIndex, main = paste("Identified peaks with SNR >", SNR.Th))
The simplest local maximum detection using a sliding window searches for maxima in a window of a given size, and slides that window across the signal, shifting it one position at a time.
localMaximum(x, winSize = 5)localMaximum(x, winSize = 5)
x |
a vector represents a signal profile |
winSize |
the slide window size, 5 by default. |
The default implementation found here shifts the window by half of its size instead of by one position at a time. This makes the implementation faster, at the expense of not being able to detect peaks that are too close to each other, if they appear in some positions with respect to the windows.
Additionally, this implementation removes all instances of peaks found at a distance less than the window size
Experimentally, we are exploring other algorithms for local maxima detection.
These algorithms can be chosen setting the "MassSpecWavelet.localMaximum.algorithm"
option. See the "Finding local maxima" vignette for further details.
Return a vector with the same length of the input x. The position of
local maximum is set as 1L, 0L else where.
Pan Du and Sergio Oller
x <- rnorm(200) lmax <- localMaximum(x, 5) maxInd <- which(lmax > 0) plot(x, type = "l") points(maxInd, x[maxInd], col = "red")x <- rnorm(200) lmax <- localMaximum(x, 5) maxInd <- which(lmax > 0) plot(x, type = "l") points(maxInd, x[maxInd], col = "red")
mexh(x)mexh(x)
x |
where to evaluate the mexican hat |
A vector of the same length as x with the corresponding values
x <- seq(-8, 8, length.out = 256) mexh(x)x <- seq(-8, 8, length.out = 256) mexh(x)
Match m/z index to m/z value with a certain error range
mzInd2vRange(mzInd, error = 0.003)mzInd2vRange(mzInd, error = 0.003)
mzInd |
a vector of m/z index |
error |
error range |
return a vector of sorted m/z values
Pan Du
Match m/z value to m/z index with a certain error range
mzV2indRange(mzV, error = 0.003)mzV2indRange(mzV, error = 0.003)
mzV |
a vector of m/z value |
error |
error range |
return a vector of sorted m/z indexes
Pan Du
This function is a wrapper of cwt(),
getLocalMaximumCWT(), getRidge(),
identifyMajorPeaks()
peakDetectionCWT( ms, scales = c(1, seq(2, 30, 2), seq(32, 64, 4)), SNR.Th = 3, nearbyPeak = TRUE, peakScaleRange = 5, amp.Th = 0.01, minNoiseLevel = amp.Th/SNR.Th, ridgeLength = 24, peakThr = NULL, tuneIn = FALSE, ..., exclude0scaleAmpThresh = FALSE, getRidgeParams = list(gapTh = 3, skip = 2) )peakDetectionCWT( ms, scales = c(1, seq(2, 30, 2), seq(32, 64, 4)), SNR.Th = 3, nearbyPeak = TRUE, peakScaleRange = 5, amp.Th = 0.01, minNoiseLevel = amp.Th/SNR.Th, ridgeLength = 24, peakThr = NULL, tuneIn = FALSE, ..., exclude0scaleAmpThresh = FALSE, getRidgeParams = list(gapTh = 3, skip = 2) )
ms |
the mass spectrometry spectrum |
scales |
Scales of CWT. See |
SNR.Th |
SNR (Signal to Noise Ratio) threshold |
nearbyPeak |
Determine whether to include the nearby small peaks of
major peaks. |
peakScaleRange |
the scale range of the peak. larger than 5 by default. |
amp.Th |
the minimum required relative amplitude of the peak (ratio to the maximum of CWT coefficients) |
minNoiseLevel |
the minimum noise level used in computing the SNR |
ridgeLength |
the minimum highest scale of the peak in 2-D CWT coefficient matrix |
peakThr |
Minimal absolute intensity (above the baseline) of peaks to
be picked. If this value is provided, then the smoothing function
|
tuneIn |
determine whether to tune in the parameter estimation of the
detected peaks. If |
... |
other parameters used by |
exclude0scaleAmpThresh |
When computing the relative |
getRidgeParams |
A list with parameters for |
majorPeakInfo |
return of |
ridgeList |
return of |
localMax |
return
of |
wCoefs |
2-D CWT coefficient
matrix, see |
Pan Du, Simon Lin
Du, P., Kibbe, W.A. and Lin, S.M. (2006) Improved peak detection in mass spectrum by incorporating continuous wavelet transform-based pattern matching, Bioinformatics, 22, 2059-2065.
cwt(), getLocalMaximumCWT(),
getRidge(), identifyMajorPeaks()
data(exampleMS) # Detect peaks with prepared wavelets: prep_wav <- prepareWavelets(length(exampleMS)) SNR.Th <- 3 peakInfo <- peakDetectionCWT(exampleMS, prep_wav, SNR.Th = SNR.Th, exclude0scaleAmpThresh=TRUE) peakIndex <- peakInfo$majorPeakInfo$peakIndex plotPeak(exampleMS, peakIndex, main = paste("Identified peaks with SNR >", SNR.Th)) SNR.Th <- 3 peakInfo <- peakDetectionCWT(exampleMS, SNR.Th = SNR.Th) majorPeakInfo <- peakInfo$majorPeakInfo peakIndex <- majorPeakInfo$peakIndex plotPeak(exampleMS, peakIndex, main = paste("Identified peaks with SNR >", SNR.Th)) ## In some cases, users may want to add peak filtering based on the absolute peak amplitude peakInfo <- peakDetectionCWT(exampleMS, SNR.Th = SNR.Th, peakThr = 500) majorPeakInfo <- peakInfo$majorPeakInfo peakIndex <- majorPeakInfo$peakIndex plotPeak(exampleMS, peakIndex, main = paste("Identified peaks with SNR >", SNR.Th))data(exampleMS) # Detect peaks with prepared wavelets: prep_wav <- prepareWavelets(length(exampleMS)) SNR.Th <- 3 peakInfo <- peakDetectionCWT(exampleMS, prep_wav, SNR.Th = SNR.Th, exclude0scaleAmpThresh=TRUE) peakIndex <- peakInfo$majorPeakInfo$peakIndex plotPeak(exampleMS, peakIndex, main = paste("Identified peaks with SNR >", SNR.Th)) SNR.Th <- 3 peakInfo <- peakDetectionCWT(exampleMS, SNR.Th = SNR.Th) majorPeakInfo <- peakInfo$majorPeakInfo peakIndex <- majorPeakInfo$peakIndex plotPeak(exampleMS, peakIndex, main = paste("Identified peaks with SNR >", SNR.Th)) ## In some cases, users may want to add peak filtering based on the absolute peak amplitude peakInfo <- peakDetectionCWT(exampleMS, SNR.Th = SNR.Th, peakThr = 500) majorPeakInfo <- peakInfo$majorPeakInfo peakIndex <- majorPeakInfo$peakIndex plotPeak(exampleMS, peakIndex, main = paste("Identified peaks with SNR >", SNR.Th))
Plot the local maximum matrix of 2-D CWT coefficients returned by
getLocalMaximumCWT()
plotLocalMax( localMax, wCoefs = NULL, range = c(1, nrow(localMax)), colorMap = "RYB", main = NULL, cex = 3, pch = ".", ... )plotLocalMax( localMax, wCoefs = NULL, range = c(1, nrow(localMax)), colorMap = "RYB", main = NULL, cex = 3, pch = ".", ... )
localMax |
local maximum matrix of 2-D CWT coefficients returned by
|
wCoefs |
2-D CWT coefficients |
range |
plot range of m/z index |
colorMap |
the colormap used in plotting the points |
main |
parameter of |
cex |
parameter of |
pch |
parameter of |
... |
other parameters of |
No value is returned; this function is called for its side effects (plot).
Pan Du
data(exampleMS) scales <- seq(1, 64, 3) wCoefs <- cwt(exampleMS[5000:11000], scales = scales, wavelet = "mexh") localMax <- getLocalMaximumCWT(wCoefs) plotLocalMax(localMax)data(exampleMS) scales <- seq(1, 64, 3) wCoefs <- cwt(exampleMS[5000:11000], scales = scales, wavelet = "mexh") localMax <- getLocalMaximumCWT(wCoefs) plotLocalMax(localMax)
Plot the identified peaks over the spectrum. The identified peaks are
returned by peakDetectionCWT() or
identifyMajorPeaks()
plotPeak( ms, peakIndex = NULL, mz = 1:length(ms), range = c(min(mz), max(mz)), method = c("p", "l"), main = NULL, log = "", ... )plotPeak( ms, peakIndex = NULL, mz = 1:length(ms), range = c(min(mz), max(mz)), method = c("p", "l"), main = NULL, log = "", ... )
ms |
the MS spectrum |
peakIndex |
m/z indexes of the identified peaks |
mz |
m/z value correspond to m/z index |
range |
the plot range of m/z value |
method |
plot method of the identified peaks. method 'p' plot circles on the peaks; method 'l' add vertical lines over the peaks. |
main |
parameter of |
log |
parameter of |
... |
other parameters of |
No value is returned; this function is called for its side effects (plot).
Pan Du
peakDetectionCWT(), identifyMajorPeaks()
data(exampleMS) SNR.Th <- 3 peakInfo <- peakDetectionCWT(exampleMS, SNR.Th = SNR.Th) majorPeakInfo <- peakInfo$majorPeakInfo peakIndex <- majorPeakInfo$peakIndex plotPeak(exampleMS, peakIndex, main = paste("Identified peaks with SNR >", SNR.Th))data(exampleMS) SNR.Th <- 3 peakInfo <- peakDetectionCWT(exampleMS, SNR.Th = SNR.Th) majorPeakInfo <- peakInfo$majorPeakInfo peakIndex <- majorPeakInfo$peakIndex plotPeak(exampleMS, peakIndex, main = paste("Identified peaks with SNR >", SNR.Th))
Plot the ridge list returned by getRidge()
plotRidgeList( ridgeList, wCoefs = NULL, range = NULL, colorMap = "RYB", main = NULL, pch = ".", cex = 3, ... )plotRidgeList( ridgeList, wCoefs = NULL, range = NULL, colorMap = "RYB", main = NULL, pch = ".", cex = 3, ... )
ridgeList |
returned by |
wCoefs |
2-D CWT coefficients |
range |
plot range of m/z index |
colorMap |
colorMap to plot the points of local maximum |
main |
parameter of |
pch |
parameter of |
cex |
parameter of |
... |
other parameters of |
No value is returned; this function is called for its side effects (plot).
Pan Du
data(exampleMS) scales <- seq(1, 64, 3) wCoefs <- cwt(exampleMS[5000:11000], scales = scales, wavelet = "mexh") localMax <- getLocalMaximumCWT(wCoefs) ridgeList <- getRidge(localMax) plotRidgeList(ridgeList)data(exampleMS) scales <- seq(1, 64, 3) wCoefs <- cwt(exampleMS[5000:11000], scales = scales, wavelet = "mexh") localMax <- getLocalMaximumCWT(wCoefs) ridgeList <- getRidge(localMax) plotRidgeList(ridgeList)
Prepare daughter wavelets for faster CWT
prepareWavelets( mslength, scales = c(1, seq(2, 30, 2), seq(32, 64, 4)), wavelet = "mexh", wavelet_xlimit = 8, wavelet_length = 1024L, extendLengthScales = TRUE )prepareWavelets( mslength, scales = c(1, seq(2, 30, 2), seq(32, 64, 4)), wavelet = "mexh", wavelet_xlimit = 8, wavelet_length = 1024L, extendLengthScales = TRUE )
mslength |
Length of the signal to transform |
scales |
a vector represents the scales at which to perform CWT. See
the Details section. Additionally, a |
wavelet |
The wavelet base, Mexican Hat by default. User can provide
wavelet |
wavelet_xlimit |
The mother wavelet will be evaluated between these limits. Ignored if |
wavelet_length |
The number of points of the mother wavelet. Ignored if |
extendLengthScales |
A logical value. If the signal is too short, we may
need to pad it to convolve it with larger daughter wavelets. Set this to |
A prepared_wavelets object.
cwt
x <- runif(2000) scales <- c(1, 2, 4, 8) prep_wavelets <- prepareWavelets(length(x), scales = scales) wCoefs <- cwt(x, prep_wavelets)x <- runif(2000) scales <- c(1, 2, 4, 8) prep_wavelets <- prepareWavelets(length(x), scales = scales) wCoefs <- cwt(x, prep_wavelets)
Smooth (denoise) the spectrum by DWT (Discrete Wavelet Transform)
smoothDWT( ms, nLevel = 6, wf = "la8", localNoiseTh = seq(1, 0, by = -0.2), localWinSize = 500, globalNoiseTh = 0.75, smoothMethod = c("soft", "hard"), method = c("dwt", "modwt") )smoothDWT( ms, nLevel = 6, wf = "la8", localNoiseTh = seq(1, 0, by = -0.2), localWinSize = 500, globalNoiseTh = 0.75, smoothMethod = c("soft", "hard"), method = c("dwt", "modwt") )
ms |
a vector representing the mass spectrum |
nLevel |
the level of DWT decomposition |
wf |
the name of wavelet for DWT |
localNoiseTh |
local noise level threshold |
localWinSize |
local window size for estimate local noise threshold |
globalNoiseTh |
global noise level threshold |
smoothMethod |
the method used for denoising. 'hard' means keeping the dwt coefficients higher than the threshold unchanged; "soft" means the dwt coefficients higher than the threshold were subtracted by the threshold. |
method |
'dwt' or 'modwt' used for decomposition |
return the smoothed mass spectrum with the 'detail' component of DWT as an attribute 'detail'.
Pan Du
Based on the identified peak position, more precise estimation of the peak
information, i.e., peak position and peak scale, can be got by this
function. The basic idea is to cut the segment of spectrum near the
identified peaks, and then do similar procedures as
peakDetectionCWT(), but with more detailed scales around the
estimated peak scale.
tuneInPeakInfo( ms, majorPeakInfo = NULL, peakIndex = NULL, peakScale = NULL, maxScale = 128, ... )tuneInPeakInfo( ms, majorPeakInfo = NULL, peakIndex = NULL, peakScale = NULL, maxScale = 128, ... )
ms |
the mass spectrometry spectrum |
majorPeakInfo |
return of |
peakIndex |
the m/z index of the identified peaks |
peakScale |
the scales of the identified peaks |
maxScale |
the maximum scale allowed for the peak |
... |
other parameters of used by |
The majorPeakInfo or peakIndex and peakScale must be provided.
peakCenterIndex |
the updated peak center m/z index |
peakScale |
the updated peak scale |
peakValue |
the corresponding peak value |
Pan Du
Du, P., Kibbe, W.A. and Lin, S.M. (2006) Improved peak detection in mass spectrum by incorporating continuous wavelet transform-based pattern matching, Bioinformatics, 22, 2059-2065.
data(exampleMS) SNR.Th <- 3 peakInfo <- peakDetectionCWT(exampleMS, SNR.Th = SNR.Th) majorPeakInfo <- peakInfo$majorPeakInfo betterPeakInfo <- tuneInPeakInfo(exampleMS, majorPeakInfo) plot(500:length(exampleMS), exampleMS[500:length(exampleMS)], type = "l", log = "x") abline(v = betterPeakInfo$peakCenterIndex, col = "red")data(exampleMS) SNR.Th <- 3 peakInfo <- peakDetectionCWT(exampleMS, SNR.Th = SNR.Th) majorPeakInfo <- peakInfo$majorPeakInfo betterPeakInfo <- tuneInPeakInfo(exampleMS, majorPeakInfo) plot(500:length(exampleMS), exampleMS[500:length(exampleMS)], type = "l", log = "x") abline(v = betterPeakInfo$peakCenterIndex, col = "red")