Issue with the .plotChromPeakDensity function?

[melpetera]

Hi,

I have noticed strange results in the past about plots obtained for groupChromPeaks regarding the maxFeatures and minFraction parameters of the PeakDensity method. As somebody reached to me to ask a question about that yesterday, I decided to properly check that and I may have spotted the issue.

The initial problem

Changing the maxFeatures parameter values does not affect the plot obtained with the plotChromPeakDensity function: always the same number of "grey rectangles". For example: 15 grey rectangles plotted on a given slice eventhough maxFeatures is set to 10.

Other issue observed: the minFraction also behave strangely on the grey rectangles -> some are plotted when the number of "points" concerned by the density peak is higher than the number of samples corresponding to the threshold implied by the minFraction parameter, although the concerned number of samples actually do not reach the threshold (as some samples have several "points"). Example: {4 samples analysed, minFraction of 0.4} -> a grey rectangle is plotted for density peaks composed of 3 points from one single sample (with no points on any of the other 3 samples).

The function used to obtain the plots

Here this issue has been observed specifically when using the plotChromPeakDensity function. To note, the tables obtained with the groupChromPeaks function are consistant with the changes of maxFeatures and minFraction values, so it is only a plot issue.

My investigation trying to understand why this occurs

I first looked at how the groupChromPeak function was plotting these grey rectangles. Arriving at the .group_peaks_density function, everything seems all right: the while loop is properly handling the two parameters for the table processing, and the plot code only highlights the areas that have been retained through the while loop. So by using the groupChromPeak function to get the plots, it should be ok.

So then I looked at the plotChromPeakDensity function (that was actually used when facing the issue), and arrived at the .plotChromPeakDensity function. As far as I understand, the "simulate" parameter is meant to enable to get plots that are consistant with the groupChromPeaks parameter one would use, by only plotting grey squares for peak groups retained by the conditions implied by the parameter choice. For that, the code is inspired by the .group_peaks_density code, in particular we find the while loop with the maxFeatures and minFraction parameters being taken into account. As the objects were not exactly the same, some adaptations have been made, and I would guess this is at that moment that the issue was introduced.

To what I see, I would say the following:

• Concerning the maxFeatures parameter, snum is defined before the while loop to count the number of features included, but the increment has not been added in the loop.
• Concerning the minFraction parameter, when comparing the testing instruction to the one from the original .group_peaks_density function, I see that the use of unique() to restrict the count to samples and not to all points has been removed.

These two observations would explain why the issue occurs. To note, I have not tested the minSamples parameter but I would guess the same issue as the one observed with minFraction should be observed, based on the things spotted in the code.

To sum up

Currently the plotChromPeakDensity function does not take into account the maxFeatures, minFraction (and potentially minSamples) parameters properly, with the issue being potentially in the underlying .plotChromPeakDensity function, with adjutments to do in the while loop inside the if (simulate) { [...] } brackets.

Potential fixes:

• adding a snum <- snum + 1 equivalent in the while loop
• adding a unique() when getting the sample groups to create the tt object

This is as far as I could reasonably investigate. Any idea?

Best regards,
Mélanie

[jorainer]

WOW, thanks @melpetera for the thorough investigation and description of the problem!

I have two suggestions:

1. would be great if you could provide a reproducible example - ideally using e.g. the data from the vignette or some other test data already used in xcms.

2. since you already spent so much time on this issue - and have now a good understanding of the code/problem - it would be cool if you could make a pull request with the suggested fix/changes to the xcms devel branch. We could then also add you as a contributor to the package - so, it would be more than just "you reporting the problem and we fixing" :)

[melpetera]

Thank you @jorainer for your feedback.

Regarding your first suggestion

Please find bellow an example using data from faahKO.

# Highlighting the problem with the plotChromPeakDensity function
# Data from faahKO is used

# 1) Setting the plot - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

library(xcms)
library(faahKO)
library(MSnbase)

# Data loading
cdfs <- dir(system.file("cdf", package = "faahKO"), full.names = TRUE,
            recursive = TRUE)[c(1, 2, 7, 8)]
pd <- data.frame(sample_name = sub(basename(cdfs), pattern = ".CDF",
                                   replacement = "", fixed = TRUE),
                 sample_group = rep("oneclass", 4),
                 stringsAsFactors = FALSE)
faahko <- readMSData(files = cdfs, pdata = new("AnnotatedDataFrame", pd), mode = "onDisk")

# findChromPeak
xchr <- findChromPeaks(faahko, param = CentWaveParam(ppm = 5, peakwidth = c(5, 20),
                                                     prefilter = c(3, 1000), mzdiff = 0.001))

# 2) Checking the behaviour of the minFraction and maxFeatures parameters with an example - - - - - - - -

# Two contrasted groupChromPeaks regarding the minFraction and maxFeatures parameters

# all: "generous" case
pdp_subs_all <- PeakDensityParam(sampleGroups = faahko$sample_group, bw = 8, binSize = 0.015,
                                 minFraction = 0.1, maxFeatures = 100)
faahko_all <- groupChromPeaks(xchr, param = pdp_subs_all)

# restricted: more selective
pdp_subs_restricted <- PeakDensityParam(sampleGroups = faahko$sample_group, bw = 8, binSize = 0.015,
                                        minFraction = 0.4, maxFeatures = 3)
faahko_restricted <- groupChromPeaks(xchr, param = pdp_subs_restricted)

# Example of peaks arround m/z 279.000

# Several ions in the 'all' case
focus_all <- which((featureDefinitions(faahko_all)$mzmed > 278.9) &
                   (featureDefinitions(faahko_all)$mzmed < 279.1))
featureDefinitions(faahko_all)[focus_all, ]

# Less ions in the 'restricted' case, in adequacy with the minFraction and maxFeatures parameters
focus_restricted <- which((featureDefinitions(faahko_restricted)$mzmed > 278) &
                          (featureDefinitions(faahko_restricted)$mzmed < 280))
featureDefinitions(faahko_restricted)[focus_restricted, ]

# 3) Checking the plots obtained with the plotChromPeakDensity function in the two configurations - - - -

par(mfrow = c(2, 1))

# 'all' case
plotChromPeakDensity(faahko_all, param = pdp_subs_all, mz = c(278.99, 279.01))
# 'restricted' case
plotChromPeakDensity(faahko_restricted, param = pdp_subs_restricted, mz = c(278.99, 279.01))

You will obtain in the end the following plot.

The first one is obtained with a minFraction of 0.1 and a generous maxFeatures of 100. The second is obtained with a minFraction of 0.4 (i.e. need of at least 2 samples this time) and a maxFeatures of 3.

From the 10 features (i.e. all features) retained in the first example, theoritically we should discard 6 features due to the minFraction change in the second example, and by changing also the maxFeatures to 3, then only keep 3.
This is actually what is observed when looking at the table results from groupChromPeaks (see part 2 in code), with no surprise.

However, we can see that concerning the plot from plotChromPeakDensity function, (i) obviously the maxFeatures is not taken into account and (ii) the minFraction is not working properly, as it effectively discard the features with only 1 "point" but keep everything alse as long as it has 2 points, even when corresponding only to the same single sample.

Regarding your second suggestion

Thank you for the suggestion, however I have to decline. I do acknowledge it would make sence and is common pratice for its convenience, but unfortunately I am very unsure I can handle the stress that making a PR on the XCMS repository would generate for me, so "me reporting" is actually the maximum I can reasonably do right now. Sorry for that and thank you in advance for your understanding.

[jorainer]

Hi, thanks for reporting back and no worries - I did not mean to put pressure on you. it's also perfectly fine to report, thanks!

[jorainer]

I had a look at your example, thanks for providing that. Just so that I understand it correctly - the main issue I see is that minFraction does only consider the total number of chromatographic peaks, but not whether they are from the same or different samples. So, in your example, even with minFraction = 0.4 a feature is defined, although (2) peaks are identified only in a single sample. Is that correct?

If that's true it seems to be a problem of the actual correspondence analysis step, not just the plotting function.

[melpetera]

I had a look at your example, thanks for providing that. Just so that I understand it correctly - the main issue I see is that minFraction does only consider the total number of chromatographic peaks, but not whether they are from the same or different samples. So, in your example, even with minFraction = 0.4 a feature is defined, although (2) peaks are identified only in a single sample. Is that correct?

Concerning the minFraction, yes this is correct

If that's true it seems to be a problem of the actual correspondence analysis step, not just the plotting function.

In fact, this problem only occurs specifically in this plotting function. In the .group_peaks_density function, line 502 the tt object is constructed to check whether the minFraction and minSamples thresholds are reached by considering actually the number of samples (which is what the parameters are about), using the trick of counting the unique sample labels. It it works well. However, in the .plotChromPeakDensity function which is noted as adapted by the very same .group_peaks_density, the corresponding tt object line 692 has "lost" the unique(), leading to counting all the chromatographic peaks instead of all the samples. This is why, I presume, it behaves the wrong way, in the plotChromPeakDensity function only.

[jorainer]

yes, thanks for the update @melpetera - I also added now a unit test to cover this, and groupChromPeaks() works as expected - only the plotChromPeakDensity(... , simulate = TRUE) has this issue. I will look into that, thanks!

[jorainer]

I can confirm that your suggestions above work, i.e., the unique() call and incrementing snum in the while. I will include these fixes in the next version. Thanks @melpetera !

[melpetera]

Cool! 🙂
Thank you @jorainer for your time, work and understanding!