###################################################
### chunk number 1: 
###################################################
library(MassSpecWavelet)


###################################################
### chunk number 2: 
###################################################
data(exampleMS)


###################################################
### chunk number 3: 
###################################################
scales <- seq(1, 64, 2)
 wCoefs <- cwt(exampleMS, scales = scales, wavelet = "mexh")


###################################################
### chunk number 4: 
###################################################
## Plot the 2-D CWT coefficients as image (It may take a while!)
xTickInterval <- 1000
plotRange <- c(5000, 11000)
image(plotRange[1]:plotRange[2], scales, wCoefs[plotRange[1]:plotRange[2],], col=terrain.colors(256), axes=FALSE, xlab='m/z index', ylab='CWT coefficient scale', main='CWT coefficients')
axis(1, at=seq(plotRange[1], plotRange[2], by=xTickInterval))
axis(2, at=c(1, seq(10, 64, by=10)))
box()


###################################################
### chunk number 5: 
###################################################
## Attach the raw spectrum as the first column
wCoefs <- cbind(as.vector(exampleMS), wCoefs)
colnames(wCoefs) <- c(0, scales)
localMax <- getLocalMaximumCWT(wCoefs)
    


###################################################
### chunk number 6: 
###################################################
plotLocalMax(localMax, wCoefs, range=plotRange)


###################################################
### chunk number 7: 
###################################################
ridgeList <- getRidge(localMax)


###################################################
### chunk number 8: 
###################################################
plotRidgeList(ridgeList,  wCoefs, range=plotRange)


###################################################
### chunk number 9: 
###################################################
SNR.Th <- 3
nearbyPeak <- TRUE
majorPeakInfo <- identifyMajorPeaks(exampleMS, ridgeList, wCoefs, SNR.Th = SNR.Th, nearbyPeak=nearbyPeak)
## Plot the identified peaks
peakIndex <- majorPeakInfo$peakIndex


###################################################
### chunk number 10: 
###################################################
plotPeak(exampleMS, peakIndex, range=plotRange, main=paste('Identified peaks with SNR >', SNR.Th)) 


###################################################
### chunk number 11: 
###################################################
data(exampleMS)
SNR.Th <- 3
nearbyPeak <- TRUE
peakInfo <- peakDetectionCWT(exampleMS, SNR.Th=SNR.Th, nearbyPeak=nearbyPeak)
majorPeakInfo = peakInfo$majorPeakInfo
peakIndex <- majorPeakInfo$peakIndex
plotRange <- c(5000, length(exampleMS))


###################################################
### chunk number 12: 
###################################################
plotPeak(exampleMS, peakIndex, range=plotRange, log='x', main=paste('Identified peaks with SNR >', SNR.Th)) 


###################################################
### chunk number 13: 
###################################################
	peakSNR <- majorPeakInfo$peakSNR
	allPeakIndex <- majorPeakInfo$allPeakIndex


###################################################
### chunk number 14: 
###################################################
	plotRange <- c(5000, 36000)
	selInd <- which(allPeakIndex >= plotRange[1] & allPeakIndex < plotRange[2])
	plot(allPeakIndex[selInd], peakSNR[selInd], type='h', xlab='m/z Index', ylab='Signal to Noise Ratio (SNR)', log='x')
	points(peakIndex, peakSNR[names(peakIndex)], type='h', col='red')
	title('Signal to Noise Ratio (SNR) of the peaks (CWT method)')


###################################################
### chunk number 15: 
###################################################
        betterPeakInfo <- tuneInPeakInfo(exampleMS, majorPeakInfo)


###################################################
### chunk number 16: 
###################################################
	plotRange <- c(5000, 11000)
        	plot(plotRange[1]:plotRange[2], exampleMS[plotRange[1]:plotRange[2]], type='l', log='x', xlab='m/z Index', ylab='Intensity')
        	abline(v=betterPeakInfo$peakCenterIndex, col='red')