Metabolomics Society Forum

Software => XCMS => R => XCMS - FAQ => Topic started by: hpbenton on August 19, 2011, 10:55:59 AM

Title: fillPeaks has zeros ?
Post by: hpbenton on August 19, 2011, 10:55:59 AM
fillPeaks() will add intensities for peaks not observed in a certain sample, but in some/most of the other samples of a smple class.
Some people have observed, that even after fillPeaks(), some values are still "0" or even "NaN" (NotANumber).
fillPeaks() is trying to (blindly) extract/integrate intensities from those raw files, where the peaks are "missing" after the peak picking step (matchedFilter or centWave).
The information, where to look, is taken from the other samples where the peak was observed, i.e. the rtmin/rtmax and mzmin/mzmax after group()ing. Two cases where this will fail are most prominent:
[attachment=1:2z4nlwmv]300px-ZeroPeaks-26.png[/attachment:2z4nlwmv] raw data just above raw data fillPeaks looking "beyond" raw data: The red EIC ends, because the peak is located at the end of the run, and the retention time has been corrected previously.

[attachment=1:2z4nlwmv]300px-ZeroPeaks-26.png[/attachment:2z4nlwmv] fillPeaks looking "beyond" raw data: The red EIC ends, because the peak is located at the end of the run, and the retention time has been corrected previously.


This is the code used:
Code: [Select]
library(xcms)
load("Anal10909_SP.Rdata")
xs <- fill.ret.Anal10909_SP

# Shortcut: aliases for groups
gr <- groups(xs)
gv <- groupval(xs, value="into")

# Shortcut: get offending group(s):
nanidx <- which(is.na(gv), arr.ind=TRUE)
zeridx <- which(gv==0, arr.ind=TRUE)

# Plot offending group(s):
e <- getEIC(xs, groupidx=nanidx[1], )
# Bad NaN guy in red:
col=rep("grey", length(sampnames(xs)))
col[nanidx[2]] <- "red"
plot(e, xs, col=col, sleep=1)
Code: [Select]
## PlotRaws
xr <- xcmsRaw(filepaths(xs)[nanidx[2]])
rtcor <- xs@rt$corrected
xr@scantime <- rtcor[[nanidx[2]]]

pdf(file="zeroPeaks.pdf")
par(mfrow=c(2,1))

for (i in 1:length(zeridx[,1])) {
    # EICs
    e <- getEIC(xs, groupidx=zeridx[i,1], )
    plot(e, xs, col=as.character(col[zeridx[i,1],]), sleep=1)

    ## PlotRaws
    xr <- xcmsRaw(filepaths(xs)[zeridx[i,2]])
    currentGroup <- gr[zeridx[i,1], ]

    ## Disable warnings for empty regions

    suppressWarnings(plotRaw(xr,
                            massrange=currentGroup[c("mzmin","mzmax")]+c(-0.1,+0.1),
                            timerange=currentGroup[c("rtmin","rtmax")]+c(-5,+5),
                            log=TRUE,
                            title=paste("Group/Samp/mz/RT: ", paste(zeridx[i,], collapse=" "), paste(round(currentGroup[c("mzmin","mzmax","rtmin","rtmax")], digits=2), collapse=" "))))

    rect(currentGroup["rtmin"],
        currentGroup["mzmin"],
        currentGroup["rtmax"],
        currentGroup["mzmax"], border="grey")
}
dev.off()

[attachment deleted by admin]