rformassspectrometry / Spectra

Low level infrastructure to handle MS spectra
https://rformassspectrometry.github.io/Spectra/
37 stars 25 forks source link

implement a filter removing Orbitrap (FFT) artefact #223

Closed jorainer closed 1 year ago

jorainer commented 2 years ago

Implement a filter to remove Orbitrap artefacts (see https://github.com/sneumann/xcms/issues/94 from @stanstrup).

sneumann commented 2 years ago

Just for posterity, here is a code snippet @stanstrup shared during the last metaRbolomics hackathon:

# Background: https://github.com/sneumann/xcms/issues/94#issuecomment-336829325

## Loading an Orbitrap spectrum with artefacts.
data(fft_spectrum)
length(mz(fft_spectrum)[[1]])
plotSpectra(fft_spectrum, xlim = c(264.5, 265.5))

plotSpectra(fft_spectrum, xlim = c(264.5, 265.5), ylim = c(0, 5e6))

fft_spectrum_filtered <- filterFourierTransformArtefacts(fft_spectrum)

fft_spectrum_filtered
length(mz(fft_spectrum_filtered)[[1]])
plotSpectra(fft_spectrum_filtered, xlim = c(264.5, 265.5), ylim = c(0, 5e6))

# Using a few examples peaks in your data you can optimize the parameters
fft_spectrum_filtered <- filterFourierTransformArtefacts(fft_spectrum,
                                                halfWindowSize = 0.2,
                                                  threshold = 0.005,
                                                  keepIsotopes = TRUE,
                                                  maxCharge = 5,
                                                  isotopeTolerance = 0.005
                                                  )

fft_spectrum_filtered
length(mz(fft_spectrum_filtered)[[1]])
plotSpectra(fft_spectrum_filtered, xlim = c(264.5, 265.5), ylim = c(0, 5e6))

# You can also do it on a whole file and write it back
# get orbitrap data
download.file("https://www.ebi.ac.uk/metabolights/ws/studies/MTBLS469/download/4cc5d820-dc5d-4766-8112-7a05f74acef4?file=AV_01_v2_male_arm1_juice.mzXML", "AV_01_v2_male_arm1_juice.mzXML")

# read data
data <- Spectra("AV_01_v2_male_arm1_juice.mzXML")

# See a problematic peak
# get example spectrum
extracted_spectrum <- data[195]

# Zoom to what should be a single peak
plotSpectra(extracted_spectrum, xlim = c(264, 266))

# Zoom intensities to clearly see the garbage should peaks
plotSpectra(extracted_spectrum, xlim = c(264, 266), ylim = c(0, 5e6))

data_filtered <- filterFourierTransformArtefacts(data,
                                                halfWindowSize = 0.2,
                                                  threshold = 0.005,
                                                  keepIsotopes = TRUE,
                                                  maxCharge = 5,
                                                  isotopeTolerance = 0.005
                                                  )
# Now fixed data
# get example spectrum
extracted_spectrum <- data_filtered[195]

# Zoom to what should be a single peak
plotSpectra(extracted_spectrum, xlim = c(264, 266))

# Zoom intensities to clearly see the garbage should peaks
plotSpectra(extracted_spectrum, xlim = c(264, 266), ylim = c(0, 5e6))

export(data_filtered, file = "new.mzML", MsBackendMzR())
jorainer commented 2 years ago

Thanks @sneumann - this is actually used as an example for the new filter function.