Inject samples and fill peaks

[stanstrup]

It can be useful to have blank samples in your final peaktable such that you can filter features that are not much higher in samples than in blanks.

But pre-processing blank samples gives major problems with retcor (at least with loess) since there are too few well-behaved peaks.

So one solution I see is if was possible to inject samples in the xcmsSet object and do a blind fillpeaks. The problem here is the retention time correction. Supposedly you'd want the inject samples to be corrected as much as possible.
What I have done so far is get the "average RT correction" and use that for the blank samples. More sophisticated idea would be to interpolate based on runorder.

If it is of any use here is my current workaround:

# Find blank files --------------------------------------------------------
files_blank <- factors %>% 
                  filter(ionmode == c("pos","neg")[settings$general$mode]) %>% 
                  filter(grepl("Blank",sample_type)) %>% 
                  extract2("filepath")


# append phenoData and filepaths ------------------------------------------
add_phenoData <- files_blank %>% 
                  basename %>% 
                  gsub(".mzML","",.) %>% 
                  {data.frame(class=rep("Blank",length(.)),row.names=.)} %>% 
                  {rbind(phenoData(xset_g_r_g_fill),.)}


add_filepaths <- files_blank %>% 
                  {c(filepaths(xset_g_r_g_fill),.)}



# Get mean rt correction for each scan ------------------------------------
# this assumes same number of scans in all files
# should interpolate to common scale

rt_diff <- xset_g_r_g_fill %>% 
        {map2(.@rt$raw,.@rt$corrected, ~..2-..1)} %>% 
        lapply(unname) %>% 
        as.data.frame %>% 
        t

cor_mean <- apply(rt_diff,2,mean)



# Get scan times for blank files ------------------------------------------
raw_blanks <- files_blank %>% 
                  readMSData(mode = "onDisk")


fData <- fData(raw_blanks)


add_rawRT <- fData %>% 
              as.tbl %>% 
              select(fileIdx, seqNum, retentionTime) %>% 
              spread(fileIdx, retentionTime) %>% 
              select(-seqNum) %>% 
              as.list

# Get correct RTs
add_corRT <- add_rawRT %>% 
             map(~..1+cor_mean)



# Add the fixed data ------------------------------------------------------
xset_g_r_g_fill_blanks <- xset_g_r_g_fill
xset_g_r_g_fill_old <- xset_g_r_g_fill 


phenoData(xset_g_r_g_fill_blanks) <- add_phenoData
filepaths(xset_g_r_g_fill_blanks) <- add_filepaths

xset_g_r_g_fill_blanks@rt$raw <- c(xset_g_r_g_fill@rt$raw, add_rawRT)
xset_g_r_g_fill_blanks@rt$corrected <- c(xset_g_r_g_fill@rt$corrected, add_rawRT)




# Do new fillPeaks --------------------------------------------------------

BPPARAM_fillpeaks <- SnowParam(min(detectCores()-2,4), progressbar = TRUE)
xset_g_r_g_fill_blanks_fill <- do.call(fillPeaks,list(object = xset_g_r_g_fill_blanks,BPPARAM = BPPARAM_fillpeaks))

xset_g_r_g_fill <- xset_g_r_g_fill_blanks_fill

[jorainer]

Oh man, I really have trouble reading the code with so many pipes ;)

Question: why do you want to have the blank samples included in the same object? I would rather want to have them separate. One could then do, as you did above, do a sort of fill peaks from (one or more) blank samples using the feature definitions from the grouped chromatographic peaks (i.e. the features). One would then get a per-feature intensity measurement in the blank.

In terms of functionality, we could expose part of the internal code for the fillChromPeaks that does nothing else than integrating signals from pre-defined mz-rt regions. The latter would come from the featureDefinitions and the function could be applied to a simple OnDiskMSnExp object containing only the blank samples.

You could not do the retention time correction on the blanks - but if it was a real blank there should ideally not be any signal, so I would think not adjusting the retention time should be OK too.

[sneumann]

Blanks are far from being blank. That's why people are interested in them.
Getting the run-order is easy if you have a Bruker instrument, as it adds the sequence
to the filename. Trouble arises if you have multiple instruments. I'd recommend
to do processing on a per-analytical LC/MS setup and integrate across instruments
afterwards. Yours, Steffen

[stanstrup]

@jotsetung Yes what you said would do the trick. I only put them in the same object because that was the easiest way I could think of to do fillPeaks.

Still if it was possible to be a little smart about the RTs it would be good. Perhaps you can extract the correction and apply that to the OnDiskMSnExp directly (either the average correction as I did or interpolate from sample order)?
Run order can also be gotten from the startTimeStamp.

As @sneumann said blanks are not really blank. Contaminants both have persistent presence in which case RT correction doesn't matter, but some also appear as peaks. Once you dig there are many of those.
I removed 1/3 of my peaks by requiring 3x as much signal in my real samples compared to blanks.

[jorainer]

OK, so let's make feature request(s) from this.

[stanstrup]

An alternative could be to allow exclusion+interpolation from sample order of specific samples in retcor. So simply assume the user have ordered the files before XCMSset.

This might be simpler/cleaner for the user and I guess also to program.