getEIC gives wrong result for zero intensity scans

[stanstrup]

XCMS version: 1.49.2

I was writing a function to make EICs for many ranges (I wanted to use the raw data and rawEIC supports only one, but my function turned out rather slow).
When I get went to getEIC (thinking I could disable the profile matrix with step = 0), that is faster, and compared results I noticed that they were not the same as my manually calculated results.

So I investigated and it looks like getEIC does something wrong sometimes with scans where there are no mz peaks in range.

I shall attempt to explains below. My test file: xraw.zip (an RDS so use readRDS)

I tried both getEICOld and getEICNew but they give same result.

Range:

a = 101.07781
b = 101.08387

Old:

opts <- options()$BioC
opts$xcms$getEIC.method <- "getEICOld"
options(BioC=opts)

EIC_xcms_old <- getEIC( xraw,
                        cbind(mzmin = a, mzmax = b),    
                        cbind(rep(0,length(a)), mzmax = rep(Inf,length(b))  , 
                              step=0)    
                        )

New:

opts <- options()$BioC
opts$xcms$getEIC.method <- "getEICNew"
options(BioC=opts)

EIC_xcms_new <- getEIC( xraw,
                        cbind(mzmin = a, mzmax = b),    
                        cbind(rep(0,length(a)), mzmax = rep(Inf,length(b))  , 
                              step=0)    
)

My own function:

EIC_manual <- get_EICs(xraw, as.data.frame(cbind(mz_lower = a, mz_upper = b)))

Now lets look at the output.
EIC_manual (empty scans a dropped, which was a bug but instructive here):

[[1]]
# A tibble: 62 x 3
    scan    scan_rt intensity
   <int>      <dbl>     <dbl>
1      1 0.04880000  1759.401
2      2 0.05708333  1391.995
3      3 0.06536667  1404.903
4      4 0.07365000  1596.921
5      6 0.09021667  1986.700
6      7 0.09850000  1135.041
7      9 0.11506667  1445.067
8     10 0.12335000  1491.382
9     11 0.13163333  1694.377
10    12 0.13991667  1590.738
# ... with 52 more rows

So scan 5 and 8 is empty.

head(EIC_xcms_old@eic$xcmsRaw[[1]], 10):

         rt intensity
 [1,] 2.928  1759.401
 [2,] 3.425  1391.995
 [3,] 3.922  1404.903
 [4,] 4.419  1596.921
 [5,] 4.916  1982.329
 [6,] 5.413  1986.700
 [7,] 5.910  1135.041
 [8,] 6.407  1577.644
 [9,] 6.904  1445.067
[10,] 7.401  1491.382

head(EIC_xcms_new@eic$xcmsRaw[[1]], 10):

         rt intensity
 [1,] 2.928  1759.401
 [2,] 3.425  1391.995
 [3,] 3.922  1404.903
 [4,] 4.419  1596.921
 [5,] 4.916  1982.329
 [6,] 5.413  1986.700
 [7,] 5.910  1135.041
 [8,] 6.407  1577.644
 [9,] 6.904  1445.067
[10,] 7.401  1491.382

See that getEIC puts an intensity also for scan 5 and 8 but otherwise agree on the intensity.
I checked who is rigt looking at the scan directly which seems to confirm that my function is right and getEIC is somehow wrong:

scandata <- getScan(xraw,8)

between <- scandata[,1] > a & scandata[,1] < b

sum(between) # --> 0
sum(scandata[between,"intensity"]) # --> 0

Any idea why this is happening? Am I misunderstanding something?

I thought this might have to do with the profile matrix but it doesn't appear so since I disabled that with step=0 both in the xcmsRaw and in getEIC.

rawEIC gives the correct result too.

[jorainer]

getEIC works on the profile matrix, not the actual raw data while rawEIC directly extracts the values from the raw data. That might be the explanation (without having looked in detail yet).

This was completely confusing to me from the beginning and in the (not so far) future I would like to clear all these things/confusions.

[stanstrup]

Yeah I understand that. But what behaviour should I then expect with step / profStep = 0? What is it doing then?

[jorainer]

Funny, with my xcms (version 1.49.3) I cannot use step = 0; that throws an error (due to mass <- seq(floor(mrange[1]/step)*step, ceiling(mrange[2]/step)*step, by = step)). So getEIC will always use the profile matrix as it would not accept a step which is 0.

I wouldn't spent too much time with that method, as this is one of the candidate-methods that have to be fixed.
My goal will be to have one method and that extracts raw EIC data. Which could then, in a second step be binned.

[jorainer]

I'm closing the issue now - feel free to re-open if needed.