In nathaniel-mahieu/centWaveP: CWT Based Chromatographic Peak Detection

library(knitr)
opts_chunk$set(fig.width=11, fig.height=12, base64_images=F, fig.align="center", message=F, warning=F, fig.path='figure_comparison/')

xcms::findPeaks.centWave vs centWaveP

library(centWaveP)
library(xcms)
library(gridExtra)
library(ggplot2)

set.seed(21)

centWaveP

file = "../../../Projects/Datasets/2nm118a_nh2_negative_ecoli_1213/2NM118A_HILICAnnoWgN14_2NM111G_12.mzXML"

xr = xcmsRaw(file, profstep=0)

roi.l = cent.xr(xr, ppm = 2, prefilter = c(4,0), maxskip = 5) # ppm is added to both sides
eic.l = lapply(sample(roi.l, 3000), roiToEic, xr, padding=100)
eic.noise.l = lapply(eic.l, estimateBaselineNoise, peakwidth = c(15, 60), minslope.peak = 20000)

peaks = lapply(seq_along(eic.noise.l), function(i) {
  wave(eic.noise.l[[i]], peakwidth = c(15, 70), valleywidth.min = 9, sensitivity = 1, smooth = T)
})

index = lapply(seq(peaks), function(i) {
  wns = peaks[[i]]$descent.fold.above.descentbaseline
  cbind(roi = rep(i, length(wns)), wsn = wns, num = seq_along(wns))
  }) %>% do.call(what = rbind)

nrow(index)

plotWavePeak = function(i, eic, peaks) {

  metadata = data.frame(peaks[i,,drop=F])

  cols = c("Local SD"="red","Wavelet Bounds"="blue","Descent Bounds"="green", "Baseline" = "orchid", "Smoothed EIC" = "orange", "EIC" = "black")

  theme_nate <- function(base_size = 16)
  {
    theme_grey(base_size = base_size) %+replace%
      theme(
        strip.background = element_rect(colour="grey95"), 
        strip.text = element_text(base_size*0.9, colour="grey40", face="plain"),

        panel.background = element_blank(),
        panel.grid.major = element_line(colour="grey95"),
        panel.grid.minor = element_line(colour="grey99"),
        plot.title = element_text(size = rel(1.2)),

        axis.title.y = element_text(vjust=1, angle=90),
        axis.title = element_text(colour="grey40"),
        axis.line = element_line(colour="grey95"),
        axis.ticks = element_line(colour="grey95"),

        legend.position="none",
        legend.key = element_blank(),
        legend.text = element_text(colour="grey40"),
        legend.title  = element_text(colour="grey40", face="bold")
      )
  }

  range = c(metadata$wavelet.start.rt-25, metadata$wavelet.end.rt+25)
  ggplot(subset(data.frame(eic), rt < range[2] & rt > range[1])) + 
    geom_line(aes(x = rt, y = i, col = "EIC")) +
    geom_line(aes(x = rt, y = baseline, col = "Baseline")) +
    geom_line(aes(x = rt, y = i.sg, colour = "Smoothed EIC"), alpha = 0.7) +
    geom_line(aes(x = rt, y = noise.local.sd, colour = "Local SD"), alpha = 0.7) +
    geom_vline(data = metadata, aes(xintercept = wavelet.start.rt+0.03, colour="Wavelet Bounds"), alpha = 0.8)+
    geom_vline(data = metadata, aes(xintercept = wavelet.end.rt+0.06, colour="Wavelet Bounds"), alpha = 0.8)+
    geom_vline(data = metadata, aes(xintercept = wavelet.location.rt+0.09, colour="Wavelet Bounds"), alpha = 0.8) +
    geom_vline(data = metadata, aes(xintercept = descent.rtcentroid, colour="Descent Bounds"), alpha = 0.8) + 
    geom_vline(data = metadata, aes(xintercept = descent.rtmin, colour="Descent Bounds"), alpha = 0.8) +
    geom_vline(data = metadata, aes(xintercept = descent.rtmax, colour="Descent Bounds"), alpha = 0.8) + 
    ggtitle(paste(metadata$sn.sg.above.min)) + 
    scale_colour_manual(name="Colours",values=cols) + 
    theme_nate() + scale_y_continuous(breaks=NULL) + scale_x_continuous(breaks=NULL) +
    coord_cartesian(ylim=c(0,max(metadata$descent.maxo)*1.1)) + ggtitle(paste(sep=" ","df:", round(metadata$descent.fold.above.descentbaseline,1), "wf:", round(metadata$wavelet.fold.above.waveletbaseline, 1)))
}

SN < 1

ps0 = lapply(sample(which(index[,"wsn"] < 1), 16, replace = T), function(x) {
  plotWavePeak(unname(index[x, "num"]), eic.noise.l[[index[x,"roi"]]], peaks[[index[x,"roi"]]])
})
do.call(grid.arrange, ps0)  

1 < SN < 2

ps1 = lapply(sample(index[,"wsn"] %>% { which(. > 1  & . < 2) }, 16, replace = T), function(x) {
  plotWavePeak(unname(index[x, "num"]), eic.noise.l[[index[x,"roi"]]], peaks[[index[x,"roi"]]])
})
do.call(grid.arrange, ps1)  

2 < SN < 3

ps2 = lapply(sample(index[,"wsn"] %>% { which(. > 2  & . < 3) }, 16, replace = T), function(x) {
  plotWavePeak(unname(index[x, "num"]), eic.noise.l[[index[x,"roi"]]], peaks[[index[x,"roi"]]])
})
do.call(grid.arrange, ps2)

3 < SN < 4

ps3 = lapply(sample(index[,"wsn"] %>% { which(. > 3  & . < 5) }, 16, replace = T), function(x) {
  plotWavePeak(unname(index[x, "num"]), eic.noise.l[[index[x,"roi"]]], peaks[[index[x,"roi"]]])
})
do.call(grid.arrange, ps3)

5 < SN

ps4 = lapply(sample(index[,"wsn"] %>% { which(. > 5) }, 16, replace=T), function(x) {
  plotWavePeak(unname(index[x, "num"]), eic.noise.l[[index[x,"roi"]]], peaks[[index[x,"roi"]]])
})
do.call(grid.arrange, ps4)

xcms::findPeaks.centWave

ROI.list = xcms:::findmzROI(xr, dev = 2*1E-6, minCentroids = 4, prefilter = c(4,0), noise = 0)

peaks.xcms = findPeaks.centWave(xr, ppm = 2, peakwidth = c(15, 60), snthresh = 0, prefilter = c(4,0), noise = 0, ROI.list = sample(ROI.list, 3000))

nrow(peaks.xcms)

head(peaks.xcms)

1 < SN < 5

 ps4 = lapply(sample(peaks.xcms[,"sn"] %>% { which(. > 1 & . < 5) }, 16, replace=T), function(i) {
   scrange = c(peaks.xcms[i,"rtmin"]-30, peaks.xcms[i,"rtmax"]+30)
   scrange[scrange < 1] = 1; scrange[scrange > max(xr@scantime)] = max(xr@scantime)
    eic = rawEIC(xr, mzrange = c(peaks.xcms[i,"mzmin"], peaks.xcms[i,"mzmax"]), rtrange = scrange)

    df = data.frame(eic)

    ggplot(df) +
      geom_line(aes(x = scan, y = intensity))
   })
do.call(grid.arrange, ps4)

5 < SN < 10

 ps4 = lapply(sample(peaks.xcms[,"sn"] %>% { which(. > 5 & . < 10) }, 16, replace=T), function(i) {
   scrange = c(peaks.xcms[i,"rtmin"]-30, peaks.xcms[i,"rtmax"]+30)
   scrange[scrange < 1] = 1; scrange[scrange > max(xr@scantime)] = max(xr@scantime)
    eic = rawEIC(xr, mzrange = c(peaks.xcms[i,"mzmin"], peaks.xcms[i,"mzmax"]), rtrange = scrange)

    df = data.frame(eic)

    ggplot(df) +
      geom_line(aes(x = scan, y = intensity))
   })
do.call(grid.arrange, ps4)

10 < SN < 100

 ps4 = lapply(sample(peaks.xcms[,"sn"] %>% { which(. > 10 & . < 100) }, 16, replace=T), function(i) {
   scrange = c(peaks.xcms[i,"rtmin"]-30, peaks.xcms[i,"rtmax"]+30)
   scrange[scrange < 1] = 1; scrange[scrange > max(xr@scantime)] = max(xr@scantime)
    eic = rawEIC(xr, mzrange = c(peaks.xcms[i,"mzmin"], peaks.xcms[i,"mzmax"]), rtrange = scrange)

    df = data.frame(eic)

    ggplot(df) +
      geom_line(aes(x = scan, y = intensity))
   })
do.call(grid.arrange, ps4)

SN > 100

 ps4 = lapply(sample(peaks.xcms[,"sn"] %>% { which(. > 100) }, 16, replace=T), function(i) {
   scrange = c(peaks.xcms[i,"rtmin"]-30, peaks.xcms[i,"rtmax"]+30)
   scrange[scrange < 1] = 1; scrange[scrange > max(xr@scantime)] = max(xr@scantime)
    eic = rawEIC(xr, mzrange = c(peaks.xcms[i,"mzmin"], peaks.xcms[i,"mzmax"]), rtrange = scrange)

    df = data.frame(eic)

    ggplot(df) +
      geom_line(aes(x = scan, y = intensity))
   })
do.call(grid.arrange, ps4)

nathaniel-mahieu/centWaveP documentation built on May 23, 2019, 12:19 p.m.