rt_idx and find peaks?#open

[deevanshi243]

Hello again :)

I'm alittle confused about the arguments like rt_idx, and dt_idx- are these the residuals in the column? What do these mean specifically and how do I know what is the optimal value for them? I read the vignette and I still couldnt really understand, would you mind providing a little more elaborate explanation?

I'm applying the pwt for peak detectoion but my code isnt doing a very good job- I'm only able to ID like 2 peaks in a sample even though I can clearly see more. So i might be setting the arguments inefficiently?

Eg. this is for a chromatogram-it starts with the scales but then just stops, shows nothing after that

 ch1 <- findPeaks(
+   one_ims2_smoothed,
+   verbose = TRUE,
+   length_in_xunits = 3,
+   peakwidth_range_xunits = c(10, 25),
+   extension_factor = 0,
+   iou_overlap_threshold = 0.2,
+   debug = FALSE
+ )
ℹ Using the following scales: 1, 10, 15, 23, 25

Super grateful! :)
Thankyou!

[deevanshi243]

Also, it says it automatically uses a reference for alignment if we dont provide one, so when I align it, my spectrum is aligned well, but my chromatogram is just one decreasing straight line (in the RIC plot). What might be doing wrong in the alignment?

[zeehio]

Hello again :)

I'm alittle confused about the arguments like rt_idx, and dt_idx- are these the residuals in the column?

I need a bit more context. What function are you using where you are passing rt_idx and dt_idx?

What do these mean specifically and how do I know what is the optimal value for them? I read the vignette and I still couldnt really understand, would you mind providing a little more elaborate explanation?

rt_ and dt_ stands for retention time and drift time respectively. _idx stands for index and typically refers to a row or column in a intensity matrix, or a position in a chromatogram or IMS spectrum.

Without a proper reprex it's hard for me to help you further here.

I'm applying the pwt for peak detectoion but my code isnt doing a very good job- I'm only able to ID like 2 peaks in a sample even though I can clearly see more. So i might be setting the arguments inefficiently?

Eg. this is for a chromatogram-it starts with the scales but then just stops, shows nothing after that

 ch1 <- findPeaks(
+   one_ims2_smoothed,
+   verbose = TRUE,
+   length_in_xunits = 3,
+   peakwidth_range_xunits = c(10, 25),
+   extension_factor = 0,
+   iou_overlap_threshold = 0.2,
+   debug = FALSE
+ )
ℹ Using the following scales: 1, 10, 15, 23, 25

Super grateful! :) Thankyou!

Could you please

saveRDS(one_ims2_smoothed, "one_ims2_smoothed.rds")

And attach the one_ims2_smoothed.rds file to this issue so I can reproduce your problem in my computer?

Thanks!

[deevanshi243]

Hey, its not letting me attach a .rds file? It says not allowed, is there some other way to do this? Would a reprex for the find peaks issue help?

[zeehio]

Put it in a zip file (if less than 25MB), otherwise upload it to Google drive, Dropbox or whatever and share a link here

[deevanshi243]

one_ims2_smoothed.zip

[deevanshi243]

New Compressed (zipped) Folder.zip

Hey, if the above doesnt work: Here are two samples from my dataset--I convert them to .csv then GCIMS objects.

[zeehio]

I've been sick these days. As soon as I get fully back on track I'll do my best to check your data

[deevanshi243]

No problem and no pressure! Take care! :)

[deevanshi243]

Hey,
I hope you are feeling better! :)
did you get a chance to review it yet?
It would help me alot!

Best,
Dee

[zeehio]

Sorry I did not have the availability to look into your issue until today.

Exploring Peak detection

Since I have made some improvements that need a bit of further testing,
I have pushed them to a specific branch.

Please install that branch with:

remotes::install_github("sipss/GCIMS@feature/issue27-peakdetection")

Let’s start by loading libraries:

library(ggplot2)
library(GCIMS)

Load and plot data

We want to load and plot our data to get a feeling of how the peaks look
like.

one_ims2_smoothed <- readRDS("one_ims2_smoothed.rds")

plot(one_ims2_smoothed)

Given the previous figure, we choose the following peak width ranges.
Peak width estimation and the whole peak detection algorithm has lots of
heuristics, so these values can’t be strictly enforced right now. This
means that you may detect peaks outside this range (and you probably
will). But still, it is a good start

peak_width_min <- 5 # seconds
peak_width_max <- 25 # seconds

Find peaks

Create sensible default CWT parameters for our data.

# If this function "does not exist" please run:
#  remotes::install_github("sipss/GCIMS@feature/issue27-peakdetection")
# to install a beta branch of the package

peakDetectionCWTParams <- defaultPeakDetectionCWTParams(
  x = rtime(one_ims2_smoothed),
  peakwidth_range_xunits = c(peak_width_min, peak_width_max) 
)
peakDetectionCWTParams

## $exclude0scaleAmpThresh
## [1] TRUE
## 
## $nearbyPeak
## [1] TRUE
## 
## $nearbyWinSize
## [1] 25
## 
## $ridgeLength
## [1] 5
## 
## $excludeBoundariesSize
## [1] 2
## 
## $SNR.Th
## [1] 3
## 
## $minNoiseLevel
## [1] 0.001
## 
## $scales
## Prepared wavelets for a spectrum of length 208, using scales: 1, 5, 7, 11, 17, 25. 
## 
## attr(,"peakwidth_range_xunits")
## [1]  5 25

Now let’s run findPeaks, plot the results and change some parameters to
see the impact:

ch1 <- findPeaks(
  one_ims2_smoothed,
  verbose = TRUE,
  length_in_xunits = 3,
  peakwidth_range_xunits = c(5, 25),
  extension_factor = 0,
  iou_overlap_threshold = 0.2,
  peakDetectionCWTParams = peakDetectionCWTParams,
  debug = TRUE
)

plot(ch1) + 
  geom_vline(
    xintercept=peaks(ch1)$apex,
    color="red"
  )

Let’s lower the signal to noise ratio threshold so we detect more peaks:

peakDetectionCWTParams$SNR.Th <- 0

ch1 <- findPeaks(
  one_ims2_smoothed,
  verbose = TRUE,
  length_in_xunits = 3,
  peakwidth_range_xunits = c(5, 25),
  extension_factor = 0,
  iou_overlap_threshold = 0.2,
  peakDetectionCWTParams = peakDetectionCWTParams,
  debug = TRUE
)

# The height at 90% of the plot, so we can place labels at a convenient place:
intensity_90range <- min(intensity(ch1)) + 0.9*(max(intensity(ch1)) - min(intensity(ch1)))
intensity_40range <- min(intensity(ch1)) + 0.4*(max(intensity(ch1)) - min(intensity(ch1)))

peak_snr <- ch1@peaks_debug_info$`1`$peakdet$majorPeakInfo$peakSNR[
  names(ch1@peaks_debug_info$`1`$peakdet$majorPeakInfo$peakIndex)
]

peak_rl <- ch1@peaks_debug_info$`1`$peakdet$majorPeakInfo$peakRidgeLengthScale[
  names(ch1@peaks_debug_info$`1`$peakdet$majorPeakInfo$peakIndex)
]

plot(ch1) + 
  geom_vline(
    xintercept=peaks(ch1)$apex,
    color="red"
  ) +
  geom_label(
    data = data.frame(
      x = ch1@peaks_debug_info$`1`$peakdet$majorPeakInfo$peakIndex,
      y = intensity_90range,
      label = paste("SNR:", format(peak_snr, digits=2))
    ),
    mapping = aes(
      x = x,
      y = y,
      label = label,
    ),
    size=3,
    color = "red",
    position = position_jitter(height=4, seed=100)
  ) +
  geom_label(
    data = data.frame(
      x = ch1@peaks_debug_info$`1`$peakdet$majorPeakInfo$peakIndex,
      y = intensity_40range,
      label = paste("RidgeLength:", format(peak_rl, digits=2))
    ),
    mapping = aes(x = x, y = y,label = label),
    size=3,
    color = "blue",
    position = position_jitter(height=4, seed=100)
  ) +
  labs(caption="Now we see too many peaks, we can start to understand how to filter them out")

Two relevant parameters to filter the peaks are the Signal to Noise
Ratio threshold SNR.Th and the ridgeLength.

peakDetectionCWTParams$SNR.Th <- 0
peakDetectionCWTParams$ridgeLength <- 7
peakDetectionCWTParams$nearbyPeak <- FALSE # Change to TRUE and you will see 
                                           # some extra peaks closer to these ones

ch1 <- findPeaks(
  one_ims2_smoothed,
  verbose = TRUE,
  length_in_xunits = 3,
  peakwidth_range_xunits = c(5, 25),
  extension_factor = 0,
  iou_overlap_threshold = 0.2,
  peakDetectionCWTParams = peakDetectionCWTParams,
  debug = TRUE
)

plot(ch1) + 
  geom_vline(
    xintercept=peaks(ch1)$apex,
    color="red"
  )

This shows a bit how I played with your data, hopefully you can reproduce this analysis and learn about the most important peak detection parameters when using wavelets.

Please ask if you have further questions.