403817 - SPE-GC/MS data analysis

Chapter 1 - Select relevant files for analysis

This is an R Markdown document. Here the developed R-script is described in detail.

Code is represented like this:

# Set working directory

setwd("~\\Adviesprojecten\\2021\\HRMS PoC\\rmarkdown\\")

When user input is necessary, it is represented like this: user input required.

Initialization

First, required packages must be loaded.

# Load required packages

library("rawrr") # required for reading spectra

library("matrixStats") # required for binning

library("spectacles") # required for SNV normalization

library("pls") # required for MSC normalization

library("signal") # required for Savitzky-Golay smoothing

library("baseline") # required for modified polynomial fitting (smoothing) and Gaussian weighting (smoothing)

library("factoextra") # For PCA

library("ggplot2") # For plotting

library("ggpubr") # For plotting

library("pheatmap") # For hierarchical clustering

library("OrgMassSpecR") # For calculating spectrum similarity

Load data

Load sample information from the CSV file generated in the previous Rmarkdown document. The user can select the files of interest here.

# Load data

samples <- read.csv(file = "\\\\nwg\\dfs\\projectdata\\P403817_001\\iQxTT2016_samplefiles.csv")



# Select all samples of September 2018, Lobith

samples <- samples[grep("1809", samples$Chromatogramm.Dateiname), ] # Select all measurements of September 2018

samples <- samples[grep("LOB", samples$Chromatogramm.Dateiname), ] # Select all Lobith measurements

Load TIC’s for all selected measurements. Make sure that the correct path is entered here.

# Create empty list for all TIC's

TIC <- vector(mode = "list", length = nrow(samples))



# Load TIC's for all pilot files

for (i in 1:nrow(samples)) {

  file <- paste0("\\\\nwg\\dfs\\projectdata\\P403817_001\\RWSdata\\", samples$filepath[i])

  temp <- rawrr::readChromatogram(rawfile = file, type = "tic")

  df <- data.frame(as.numeric(temp$times), temp$intensities)

  colnames(df) <- c("time", "intensity")

  TIC[[i]] <- df

  names(TIC)[i] <- samples$Chromatogramm.Dateiname[i]

  print(i)

  rm(file, temp, df)

}



# Store list in an RDS object

saveRDS(object = TIC, file = "TIC_validation_1month_phenol.RDS")

Load MS spectra for all selected measurements. Make sure that the correct path is entered here.

# Create empty list for MS spectra

spectra <- vector(mode = "list", length = nrow(samples))

names(spectra) <- samples$Chromatogramm.Dateiname



# Read MS spectra for samples

for (k in 1:nrow(samples)) {

  file <- paste0("\\\\nwg\\dfs\\projectdata\\P403817_001\\RWSdata\\", samples$filepath[k])

  filename <- samples$Chromatogramm.Dateiname[k]

  n <- rawrr::readFileHeader(rawfile = file)$`Number of scans`

  temp <- rawrr::readSpectrum(rawfile = file, scan = c(1:n))

  spectra[[filename]] <- temp

  rm(temp, file, filename, n)

  print(k)

}



# Store list in an RDS object

saveRDS(object = spectra, file = "MSspectra_validation_1month_phenol.RDS")

Load reference data. User input required in the section below, he/she must provide the reference chromatogram and the file with reference MS spectra of the 8 internal standards.

# Read file with reference MS spectra of the 8 internal standards

IS <- readRDS(file = "IS_MSspectra_pilot.RDS")



# To be filled in by the operator, reference chromatogram used for KRetI alignment and other preprocessing steps (such as MSC)

ref.chromatogram <- "180901_LOB_06"

Define required functions

Define function to retrieve the retention times from internal standards quickly. This function is applied in Chapter 8 - PCA.

getRt <- function(all.ref.IS, chromatogram) {

  chrom.rt.IS <- data.frame(matrix(ncol = 8, nrow = 1))

  colnames(chrom.rt.IS) <- c(

    "toluene-d8", "chlorobenzene-d5", "dichlorobenzene-d4", "naphthalene-d8",

    "dibromobenzene-d4", "terbuthylazine-d5", "phenanthrene-d10", "chrysene-d12"

  )

  chrom.rt.IS$`toluene-d8` <- all.ref.IS[["toluene-d8"]][["Rt.IS"]][which(all.ref.IS[["toluene-d8"]][["Rt.IS"]]$FileName == chromatogram), "rtinmin"]

  chrom.rt.IS$`chlorobenzene-d5` <- all.ref.IS[["chlorobenzene-d5"]][["Rt.IS"]][which(all.ref.IS[["chlorobenzene-d5"]][["Rt.IS"]]$FileName == chromatogram), "rtinmin"]

  chrom.rt.IS$`dichlorobenzene-d4` <- all.ref.IS[["dichlorobenzene-d4"]][["Rt.IS"]][which(all.ref.IS[["dichlorobenzene-d4"]][["Rt.IS"]]$FileName == chromatogram), "rtinmin"]

  chrom.rt.IS$`naphthalene-d8` <- all.ref.IS[["naphthalene-d8"]][["Rt.IS"]][which(all.ref.IS[["naphthalene-d8"]][["Rt.IS"]]$FileName == chromatogram), "rtinmin"]

  chrom.rt.IS$`dibromobenzene-d4` <- all.ref.IS[["dibromobenzene-d4"]][["Rt.IS"]][which(all.ref.IS[["dibromobenzene-d4"]][["Rt.IS"]]$FileName == chromatogram), "rtinmin"]

  chrom.rt.IS$`terbuthylazine-d5` <- all.ref.IS[["terbuthylazine-d5"]][["Rt.IS"]][which(all.ref.IS[["terbuthylazine-d5"]][["Rt.IS"]]$FileName == chromatogram), "rtinmin"]

  chrom.rt.IS$`phenanthrene-d10` <- all.ref.IS[["phenanthrene-d10"]][["Rt.IS"]][which(all.ref.IS[["phenanthrene-d10"]][["Rt.IS"]]$FileName == chromatogram), "rtinmin"]

  chrom.rt.IS$`chrysene-d12` <- all.ref.IS[["chrysene-d12"]][["Rt.IS"]][which(all.ref.IS[["chrysene-d12"]][["Rt.IS"]]$FileName == chromatogram), "rtinmin"]

  return(chrom.rt.IS)

}

Define function to add classification (in this case of location). This function is applied in Chapter 2.

getClassification <- function(name) {

  if (grepl("LOB", name)) {

    a <- "LOB"

  } else if (grepl("BIM", name)) {

    a <- "BIM"

  } else if (grepl("REE", name)) {

    a <- "REE"

  } else if (grepl("WSL", name)) {

    a <- "WSL"

  } else if (grepl("GWH", name)) {

    a <- "GWH"

  } else if (grepl("EIJ", name)) {

    a <- "EIJ"

  } else {

    a <- "unknown"

  }

  return(a)

}

Define function to add classification (in this case of year). This function is applied in Chapter 8 - PCA.

getYear <- function(name) {

  if (substr(name, start = 1, stop = 2) == "19") {

    a <- "2019"

  } else if (substr(name, start = 1, stop = 2) == "20") {

    a <- "2020"

  } else if (substr(name, start = 1, stop = 2) == "21") {

    a <- "2021"

  } else if (substr(name, start = 1, stop = 2) == "18") {

    a <- "2018"

  } else if (substr(name, start = 1, stop = 2) == "17") {

    a <- "2017"

  } else if (substr(name, start = 1, stop = 2) == "16") {

    a <- "2016"

  } else {

    a <- "unknown"

  }

  return(a)

}

Chapter 2 - Localization of internal standards

Retrieve per internal standard (IS) the reference fragments (quantifier and qualifiers) used to find the retention time. In a later stage, this part should be done by the analyst, who provides the m/z of the Q1-Q4 and the script will calculate their relative intensities. The algorithm relies on the specific MS spectrum of the internal standard. In particular, a reference MS-spectrum and expected Rt of the internal standard are defined by the operator. Subsequently, the algorithm defines an Rt window, based on the Rt provided by the user, and uses a spectral similarity function to find which of the recorded spectra in the window matches the one defined by the user. User input required in the sections below, he/she must provide the retention time window where to look for the internal standard.

# Manually define Rt window (in min) to look for the IS. Tested also 2min but there seems to be no difference in performances.

Rt.window <- 3



# Initialise a list which will contain all info about the reference IS

all.ref.IS <- list()

#' Define a filtering function to remove intensities below 5% of the max

low_int <- function(MSspec) {

  MSspec <- MSspec[which(MSspec$intensity > max(MSspec$intensity, na.rm = T) * 0.05), ]

}



#' Define a function to normalize the intensities

norm_int <- function(MSintensity) {

  MSintensity <- (MSintensity / max(MSintensity, na.rm = T)) * 100

}

Prepare reference m/z of each internal standard.
Toluene-d8

# Retrieve mZ & intensity for the IS

ref.IS.spectrum <- data.frame(mZ = IS[["toluene-d8"]]$mZ, intensity = IS[["toluene-d8"]]$intensity)

ref.IS.spectrum <- low_int(ref.IS.spectrum) # remove m/z which intensity is less than 5% of the max

ref.IS.spectrum$intensity <- norm_int(ref.IS.spectrum$intensity) # normalize intensities

ref.IS.spectrum.top4 <- ref.IS.spectrum[order(ref.IS.spectrum$intensity, decreasing = T), ] # order by decreasing intensity



# Manually select most abundant and diagnostic fragments for IS 1

ref.IS.frag <- ref.IS.spectrum

ref.IS.frag$round.mz <- round(ref.IS.frag$mZ / 0.5) * 0.5

ref.IS.frag$rel.int <- (ref.IS.frag$intensity / max(ref.IS.frag$intensity)) * 100



# Manually define the Rt (in min)

expected.Rt <- 4.7



# Calculate the range

expected.Rt.range <- c(expected.Rt - Rt.window, expected.Rt + Rt.window) # define the Rt range to look for the IS

expected.Rt.range <- expected.Rt.range * 60 # transform in seconds



all.ref.IS[["toluene-d8"]]$ref.IS.frag <- ref.IS.frag

all.ref.IS[["toluene-d8"]]$RtRange <- expected.Rt.range

all.ref.IS[["toluene-d8"]]$expected.Rt <- expected.Rt

For Toluene-d8, this results in the following list:

##           mZ  intensity round.mz    rel.int

## 10  42.06339  12.929432       42  12.929432

## 11  43.06023   7.532936       43   7.532936

## 20  52.04442   5.324995       52   5.324995

## 22  54.06854   7.325357       54   7.325357

## 34  66.06665   6.231184       66   6.231184

## 38  70.07024  13.921788       70  13.921788

## 59  91.02415   5.520571       91   5.520571

## 65  97.09407   5.484367       97   5.484367

## 66  98.09047 100.000000       98 100.000000

## 67  99.12589  11.537655       99  11.537655

## 68 100.10323  60.227808      100  60.227808

## [1] 102 462

## [1] 4.7

Now, the same is done for all other internal standards.

Chlorobenzene-d5

# Retrieve mZ & intensity for the IS

ref.IS.spectrum <- data.frame(mZ = IS[["chlorobenzene-d5"]]$mZ, intensity = IS[["chlorobenzene-d5"]]$intensity)

ref.IS.spectrum <- low_int(ref.IS.spectrum) # remove mz which intensity is less than 5% of the max

ref.IS.spectrum$intensity <- norm_int(ref.IS.spectrum$intensity) # normalize intensities

ref.IS.spectrum.top4 <- ref.IS.spectrum[order(ref.IS.spectrum$intensity, decreasing = T), ] # order by decreasing intensity



# Manually select most abundant and diagnostic fragments for IS 2

ref.IS.frag <- ref.IS.spectrum

ref.IS.frag$round.mz <- round(ref.IS.frag$mZ / 0.5) * 0.5

ref.IS.frag$rel.int <- (ref.IS.frag$intensity / max(ref.IS.frag$intensity)) * 100



# Manually define the Rt (in min)

expected.Rt <- 4.97



# Calculate the range

expected.Rt.range <- c(expected.Rt - Rt.window, expected.Rt + Rt.window) # define the Rt range to look for the IS

expected.Rt.range <- expected.Rt.range * 60 # transform in seconds



all.ref.IS[["chlorobenzene-d5"]]$ref.IS.frag <- ref.IS.frag

all.ref.IS[["chlorobenzene-d5"]]$RtRange <- expected.Rt.range

all.ref.IS[["chlorobenzene-d5"]]$expected.Rt <- expected.Rt

1,4-Dichlorobenzene-d4

# Retrieve mZ & intensity for the IS

ref.IS.spectrum <- data.frame(mZ = IS[["dichlorobenzene-d4"]]$mZ, intensity = IS[["dichlorobenzene-d4"]]$intensity)

ref.IS.spectrum <- low_int(ref.IS.spectrum) # remove mz which intensity is less than 5% of the max

ref.IS.spectrum$intensity <- norm_int(ref.IS.spectrum$intensity) # normalize intensities

ref.IS.spectrum.top4 <- ref.IS.spectrum[order(ref.IS.spectrum$intensity, decreasing = T), ] # order by decreasing intensity



# Manually select most abundant and diagnostic fragments for IS 3

ref.IS.frag <- ref.IS.spectrum

ref.IS.frag$round.mz <- round(ref.IS.frag$mZ / 0.5) * 0.5

ref.IS.frag$rel.int <- (ref.IS.frag$intensity / max(ref.IS.frag$intensity)) * 100



# Manually define the Rt (in min)

expected.Rt <- 5.99



# Calculate the range

expected.Rt.range <- c(expected.Rt - Rt.window, expected.Rt + Rt.window) # define the Rt range to look for the IS

expected.Rt.range <- expected.Rt.range * 60 # transform in seconds



all.ref.IS[["dichlorobenzene-d4"]]$ref.IS.frag <- ref.IS.frag

all.ref.IS[["dichlorobenzene-d4"]]$RtRange <- expected.Rt.range

all.ref.IS[["dichlorobenzene-d4"]]$expected.Rt <- expected.Rt

Naphthalene-d8

# Retrieve mZ & intensity for the IS

ref.IS.spectrum <- data.frame(mZ = IS[["naphthalene-d8"]]$mZ, intensity = IS[["naphthalene-d8"]]$intensity)

ref.IS.spectrum <- low_int(ref.IS.spectrum) # remove mz which intensity is less than 5% of the max

ref.IS.spectrum$intensity <- norm_int(ref.IS.spectrum$intensity) # normalize intensities

ref.IS.spectrum.top4 <- ref.IS.spectrum[order(ref.IS.spectrum$intensity, decreasing = T), ] # order by decreasing intensity



# Manually select most abundant and diagnostic fragments for IS 4

ref.IS.frag <- ref.IS.spectrum

ref.IS.frag$round.mz <- round(ref.IS.frag$mZ / 0.5) * 0.5

ref.IS.frag$rel.int <- (ref.IS.frag$intensity / max(ref.IS.frag$intensity)) * 100



# Manually define the Rt (in min)

expected.Rt <- 6.8



# Calculate the range

expected.Rt.range <- c(expected.Rt - Rt.window, expected.Rt + Rt.window) # define the Rt range to look for the IS

expected.Rt.range <- expected.Rt.range * 60 # transform in seconds



all.ref.IS[["naphthalene-d8"]]$ref.IS.frag <- ref.IS.frag

all.ref.IS[["naphthalene-d8"]]$RtRange <- expected.Rt.range

all.ref.IS[["naphthalene-d8"]]$expected.Rt <- expected.Rt

1,4-Dibromobenzene-d4

# Retrieve mZ & intensity for the IS

ref.IS.spectrum <- data.frame(mZ = IS[["dibromobenzene-d4"]]$mZ, intensity = IS[["dibromobenzene-d4"]]$intensity)

ref.IS.spectrum <- low_int(ref.IS.spectrum) # remove mz which intensity is less than 5% of the max

ref.IS.spectrum$intensity <- norm_int(ref.IS.spectrum$intensity) # normalize intensities

ref.IS.spectrum.top4 <- ref.IS.spectrum[order(ref.IS.spectrum$intensity, decreasing = T), ] # order by decreasing intensity



# Manually select most abundant and diagnostic fragments for IS 5

ref.IS.frag <- ref.IS.spectrum

ref.IS.frag$round.mz <- round(ref.IS.frag$mZ / 0.5) * 0.5

ref.IS.frag$rel.int <- (ref.IS.frag$intensity / max(ref.IS.frag$intensity)) * 100



# Manually define the Rt (in min)

expected.Rt <- 6.83



# Calculate the range

expected.Rt.range <- c(expected.Rt - Rt.window, expected.Rt + Rt.window) # define the Rt range to look for the IS

expected.Rt.range <- expected.Rt.range * 60 # transform in seconds



all.ref.IS[["dibromobenzene-d4"]]$ref.IS.frag <- ref.IS.frag

all.ref.IS[["dibromobenzene-d4"]]$RtRange <- expected.Rt.range

all.ref.IS[["dibromobenzene-d4"]]$expected.Rt <- expected.Rt

Terbuthylazine-d5

# Retrieve mZ & intensity for the IS

ref.IS.spectrum <- data.frame(mZ = IS[["terbuthylazine-d5"]]$mZ, intensity = IS[["terbuthylazine-d5"]]$intensity)

ref.IS.spectrum <- low_int(ref.IS.spectrum) # remove mz which intensity is less than 5% of the max

ref.IS.spectrum$intensity <- norm_int(ref.IS.spectrum$intensity) # normalize intensities

ref.IS.spectrum.top4 <- ref.IS.spectrum[order(ref.IS.spectrum$intensity, decreasing = T), ] # order by decreasing intensity



# Manually select most abundant and diagnostic fragments for IS 6

ref.IS.frag <- ref.IS.spectrum

ref.IS.frag$round.mz <- round(ref.IS.frag$mZ / 0.5) * 0.5

ref.IS.frag$rel.int <- (ref.IS.frag$intensity / max(ref.IS.frag$intensity)) * 100



# Manually define the Rt (in min)

expected.Rt <- 8.7



# Calculate the range

expected.Rt.range <- c(expected.Rt - Rt.window, expected.Rt + Rt.window) # define the Rt range to look for the IS

expected.Rt.range <- expected.Rt.range * 60 # transform in seconds



all.ref.IS[["terbuthylazine-d5"]]$ref.IS.frag <- ref.IS.frag

all.ref.IS[["terbuthylazine-d5"]]$RtRange <- expected.Rt.range

all.ref.IS[["terbuthylazine-d5"]]$expected.Rt <- expected.Rt

Phenanthrene-d10

# Retrieve mZ & intensity for the IS

ref.IS.spectrum <- data.frame(mZ = IS[["phenanthrene-d10"]]$mZ, intensity = IS[["phenanthrene-d10"]]$intensity)

ref.IS.spectrum <- low_int(ref.IS.spectrum) # remove mz which intensity is less than 5% of the max

ref.IS.spectrum$intensity <- norm_int(ref.IS.spectrum$intensity) # normalize intensities

ref.IS.spectrum.top4 <- ref.IS.spectrum[order(ref.IS.spectrum$intensity, decreasing = T), ] # order by decreasing intensity



# Manually select most abundant and diagnostic fragments for IS 7

ref.IS.frag <- ref.IS.spectrum

ref.IS.frag$round.mz <- round(ref.IS.frag$mZ / 0.5) * 0.5

ref.IS.frag$rel.int <- (ref.IS.frag$intensity / max(ref.IS.frag$intensity)) * 100



# Manually define the Rt (in min)

expected.Rt <- 8.87



# Calculate the range

expected.Rt.range <- c(expected.Rt - Rt.window, expected.Rt + Rt.window) # define the Rt range to look for the IS

expected.Rt.range <- expected.Rt.range * 60 # transform in seconds



all.ref.IS[["phenanthrene-d10"]]$ref.IS.frag <- ref.IS.frag

all.ref.IS[["phenanthrene-d10"]]$RtRange <- expected.Rt.range

all.ref.IS[["phenanthrene-d10"]]$expected.Rt <- expected.Rt

Chrysene-d12

# Retrieve mZ & intensity for the IS

ref.IS.spectrum <- data.frame(mZ = IS[["chrysene-d12"]]$mZ, intensity = IS[["chrysene-d12"]]$intensity)

ref.IS.spectrum <- low_int(ref.IS.spectrum) # remove mz which intensity is less than 5% of the max

ref.IS.spectrum$intensity <- norm_int(ref.IS.spectrum$intensity) # normalize intensities

ref.IS.spectrum.top4 <- ref.IS.spectrum[order(ref.IS.spectrum$intensity, decreasing = T), ] # order by decreasing intensity



# Manually select most abundant and diagnostic fragments for IS 8

ref.IS.frag <- ref.IS.spectrum

ref.IS.frag$round.mz <- round(ref.IS.frag$mZ / 0.5) * 0.5

ref.IS.frag$rel.int <- (ref.IS.frag$intensity / max(ref.IS.frag$intensity)) * 100



# Manually define the Rt (in min)

expected.Rt <- 10.78



# Calculate the range

expected.Rt.range <- c(expected.Rt - Rt.window, expected.Rt + Rt.window) # define the Rt range to look for the IS

expected.Rt.range <- expected.Rt.range * 60 # transform in seconds



all.ref.IS[["chrysene-d12"]]$ref.IS.frag <- ref.IS.frag

all.ref.IS[["chrysene-d12"]]$RtRange <- expected.Rt.range

all.ref.IS[["chrysene-d12"]]$expected.Rt <- expected.Rt

Next, search in scans for reference fragments to determine the retention time of the IS in the samples. User input required in the section below, he/she can adjust the tolerance parameters.

Set some parameters

# Set the m/z tolerance to search for the quantifier in the spectra

mz.tolerance <- 0.5



# Set the tolerance for deviation from relative intensity

rel.int.tolerance <- 20



# Set the tolerance for the spectrum similarity score

score.tolerance <- 0.7

Apply loop to all chromatograms

for (k in 1:length(all.ref.IS)) {



  # Initialize an empty vector which will contain the name of each chromatogram, IS and corresponding Rt

  all.Rt.IS <- data.frame()



  # Loop over spectra, a large list containing all the chromatograms that need to be searched

  for (n in 1:length(spectra)) {



    # Initialize a vector which will contain scans which have the relevant fragment ion (quantifier)

    saved.scan.num <- NULL



    # Initialise a vector which will contain the SpectrumSimilarity score of the relevant fragment ion

    saved.scan.score <- NULL



    # Retrieve all scans and their corresponding retention times in seconds.

    # This will be used for the Rt window to look for the IS

    all.Rt.scans <- NULL

    for (y in 1:length(spectra[[n]])) {

      Rt.scans <- data.frame(rtinseconds = spectra[[n]][[y]]$rtinseconds, scan = spectra[[n]][[y]]$scan)

      all.Rt.scans <- rbind(all.Rt.scans, Rt.scans) # append the scans and Rt to a vector

    }



    # Retrieve the Rt range as given by the operator

    Range <- all.ref.IS[[k]]$RtRange



    # We use a +- 2min window, hence it can happen that the range will be below 0.

    # Below we put it to 0 should the min of the range be negative

    if (min(Range) < 0) {

      Range <- c(0, max(Range))

    }



    # Similarly we check if the max range is out of bounds

    if (max(Range) > max(all.Rt.scans$rtinseconds)) {

      Range <- c(Range[1], max(all.Rt.scans$rtinseconds))

    }



    # Retrieve the scans which correspond to the relevant Rt

    selected.Rt.scans <- all.Rt.scans[which(all.Rt.scans$rtinseconds >= Range[1] &

      all.Rt.scans$rtinseconds <= Range[2]), ]



    for (i in min(selected.Rt.scans):max(selected.Rt.scans)) {

      # retrieve the first scan of the first chromatogram to be searched for

      scan.i <- spectra[[n]][[i]]



      # retrieve the m/z of the first scan and round them

      scan.i.mz <- data.frame(mZ = scan.i$mZ, intensity = scan.i$intensity)

      scan.i.mz <- low_int(scan.i.mz) # remove intensities below 5% of the max

      scan.i.mz$intensity <- norm_int(scan.i.mz$intensity) # normalize intensities



      # Calculate spectral similarity

      require(OrgMassSpecR)



      score <- SpectrumSimilarity(cbind(scan.i.mz$mZ, scan.i.mz$intensity),

        cbind(all.ref.IS[[k]]$ref.IS.frag$mZ, all.ref.IS[[k]]$ref.IS.frag$intensity),

        t = 0.1, print.graphic = F

      )



      # Save all scores for all scans

      # if the similarity between the found spectrum and the reference IS spectra is above the threshold

      if (is.na(score)) {}

      else if (score >= score.tolerance) {

        saved.scan.num <- rbind(saved.scan.num, c(scan.i$scan, "Score equal to or above Tolerance"))

      }

      else if (score < score.tolerance & score >= 0.6) {

        saved.scan.num <- rbind(saved.scan.num, c(scan.i$scan, "Score equal to or above 0.6"))

      }

      else {}

    }



    #### Below is to handle instances where the IS cannot be found. Here we introduce two options.####

    # 1) We expand the research of the IS to all scans to see if we can find it.

    # 2) If we cannot find, then we replace everything with NA



    # Option 1)

    if (is.null(saved.scan.num)) {



      # Search through all scans

      for (i in 1:length(spectra[[n]])) {

        # retrieve the first scan of the first chromatogram to be searched for

        scan.i <- spectra[[n]][[i]]



        # retrieve the m/z of the first scan and round them

        scan.i.mz <- data.frame(mZ = scan.i$mZ, intensity = scan.i$intensity)

        scan.i.mz <- low_int(scan.i.mz) # remove intensities below 5% of the max

        scan.i.mz$intensity <- norm_int(scan.i.mz$intensity) # normalize intensities



        # Calculate spectral similarity

        require(OrgMassSpecR)



        score <- SpectrumSimilarity(cbind(scan.i.mz$mZ, scan.i.mz$intensity),

          cbind(all.ref.IS[[k]]$ref.IS.frag$mZ, all.ref.IS[[k]]$ref.IS.frag$intensity),

          t = 0.1, print.graphic = F

        )



        # Save all scores for all scans

        # if the similarity between the found spectrum and the reference IS spectra is above the threshold

        if (is.na(score)) {}

        else if (score >= score.tolerance) {

          saved.scan.num <- rbind(saved.scan.num, c(

            scan.i$scan,

            "Score >= Tolerance (All scans have been searched!)"

          ))

        }

        else if (score < score.tolerance & score >= 0.6) {

          saved.scan.num <- rbind(saved.scan.num, c(

            scan.i$scan,

            "Score >= 0.6 (All scans have been searched!)"

          ))

        }

        else {}

      }

    }



    # 2) if the saved.scan.num is still empty, then we replace it with NAs

    if (is.null(saved.scan.num)) {

      rt.is <- data.frame(FileName = names(TIC[n]), rtinseconds = NA, rtinmin = NA, scan = NA, alert = NA)

    }

    # End of options 1 and 2#



    else {

      # Transform to data.frame and make scan a number

      saved.scan.num <- data.frame(scan = saved.scan.num[, 1], alert = saved.scan.num[, 2])

      saved.scan.num$scan <- as.numeric(as.character(saved.scan.num$scan))



      # Select the scans which match the criteria

      sel.scan <- spectra[[n]][saved.scan.num$scan]



      # Initialise a vector which will contain the intensities of the mz of the quantifier (main fragment)

      # in the selected spectra

      sel.scan.int <- NULL

      for (m in 1:length(sel.scan)) {

        # create a data frame only with the relevant information

        sel.scan.m.mz <- data.frame(mZ = sel.scan[[m]]$mZ, intensity = sel.scan[[m]]$intensity)



        # find the intensity of the mz which matches the mz of the quantifier used as reference.

        # The mz needs to be +- within the defined tolerance. Once its been found, append it to

        # the previously initialised vector with all the intensities

        sel.scan.int <- rbind(

          sel.scan.int,

          max(sel.scan.m.mz[which(round(sel.scan.m.mz$mZ, 1) > round(all.ref.IS[[k]]$ref.IS.frag$mZ[1], 1) -

            mz.tolerance & round(sel.scan.m.mz$mZ, 1) <

            round(all.ref.IS[[k]]$ref.IS.frag$mZ[1], 1) + mz.tolerance), ]$intensity)

        )

      }



      # find the scan number of the spectra with the highest intensity.

      sel.scan.int.max <- sel.scan[[which.max(sel.scan.int)]]$scan



      # Use the retrieved scan number to find the Rt of the IS in question

      rt.is <- spectra[[n]][sel.scan.int.max]



      # create a data.frame with all relevant information

      rt.is <- data.frame(

        FileName = names(TIC[n]), rtinseconds = rt.is[[1]]$rtinseconds, rtinmin = rt.is[[1]]$rtinseconds / 60,

        scan = sel.scan.int.max, alert = saved.scan.num$alert[which.max(sel.scan.int)]

      )

    }



    # append data.frame to data.frame initialised at the beginning of the for loop

    all.Rt.IS <- rbind(all.Rt.IS, rt.is)

  }



  # All Rt of the IS are saved in this list, including corresponding scan number and sample name.

  all.ref.IS[[k]]$Rt.IS <- all.Rt.IS



  # Print number of iteration so that the user knows how far the calculation is

  print(k)

}



# Save list to an RDS file

saveRDS(all.ref.IS, file = "Rt_IS_validation.RDS")

For Toluene-d8, this results in the following data frame with filenames, RT’s and scannumbers:

##        FileName rtinseconds  rtinmin scan                      alert

## 1 180901_LOB_06     252.700 4.211667   64 Score â<U+0089>¥ Tolerance

## 2 180901_LOB_18     252.098 4.201633   61 Score â<U+0089>¥ Tolerance

## 3 180902_LOB_06     253.100 4.218333   66 Score â<U+0089>¥ Tolerance

## 4 180902_LOB_18     253.504 4.225067   68 Score â<U+0089>¥ Tolerance

## 5 180903_LOB_06     253.718 4.228633   69 Score â<U+0089>¥ Tolerance

## 6 180903_LOB_18     252.904 4.215067   65 Score â<U+0089>¥ Tolerance

Chapter 3 - Retention time alignment

To correct for shifts in retention time (Rt), an approach based on the calculation of retention time indexes was implemented. The retention times of the internal standards are used to align the TIC’s. First, the Rt’s of the reference chromatogram are retrieved and an empty list to store the data in is prepared.

# Get Rt's of reference TIC

ref <- getRt(all.ref.IS = all.ref.IS, chromatogram = ref.chromatogram)



# Create empty list for aligned TIC's

TIC_align <- vector(mode = "list", length = length(TIC))

names(TIC_align) <- names(TIC)

Next, the retention indexes are calculated and stored in the list.

for (z in 1:length(TIC)) {

  clipb <- TIC[[z]] # select TIC

  is <- getRt(all.ref.IS = all.ref.IS, chromatogram = names(TIC[z])) # get the retention times of the IS

  is$`toluene-d8` <- NULL # This IS is not used for KRetI alignment

  is$`terbuthylazine-d5` <- NULL # This IS is not used for KRetI alignment

  is$`chrysene-d12` <- NULL # This IS is not used for KRetI alignment



  if (!any(is.na(unlist(is)))) { # Check if all IS have been found by checking if no NA is present

    clipb$time_corrected <- NA # define new column where corrected retention time can be stored

    for (j in 1:nrow(clipb)) {

      rtX <- clipb$time[j]

      if (rtX <= ref$`dichlorobenzene-d4`) {

        clipb$time_corrected[j] <- ref$`chlorobenzene-d5` +

          ((ref$`dichlorobenzene-d4` - ref$`chlorobenzene-d5`) /

            (is$`dichlorobenzene-d4` - is$`chlorobenzene-d5`)) * (rtX - is$`chlorobenzene-d5`)

      } else if (rtX > ref$`chlorobenzene-d5` && rtX <= ref$`naphthalene-d8`) {

        clipb$time_corrected[j] <- ref$`naphthalene-d8` +

          ((ref$`chlorobenzene-d5` - ref$`naphthalene-d8`) /

            (is$`chlorobenzene-d5` - is$`naphthalene-d8`)) * (rtX - is$`naphthalene-d8`)

      } else if (rtX > ref$`naphthalene-d8` && rtX <= ref$`dibromobenzene-d4`) {

        clipb$time_corrected[j] <- ref$`dibromobenzene-d4` +

          ((ref$`naphthalene-d8` - ref$`dibromobenzene-d4`) /

            (is$`naphthalene-d8` - is$`dibromobenzene-d4`)) * (rtX - is$`dibromobenzene-d4`)

      } else if (rtX > ref$`dibromobenzene-d4`) {

        clipb$time_corrected[j] <- ref$`phenanthrene-d10` +

          ((ref$`dibromobenzene-d4` - ref$`phenanthrene-d10`) /

            (is$`dibromobenzene-d4` - is$`phenanthrene-d10`)) * (rtX - is$`phenanthrene-d10`)

      }

    }



    TIC_align[[names(TIC[z])]] <- clipb[, c("time", "time_corrected", "intensity")]

    rm(clipb, is, j, rtX)

  } else { # if NA or Inf is present perform this

    TIC_align[[names(TIC[z])]] <- NULL

  } # close if else statement



  print(z) # print number of iteration to track progress

} # close for loop



# Save object

saveRDS(object = TIC_align, file = "TIC_validation_aligned_KRetI.RDS")

rm(TIC, z) # remove variables that we don't need anymore

Chapter 4 - Binning

The TIC’s are binned using a bin size of 0.01 minutes.

delta <- 0.01 # Define bin size here

TIC_binned <- vector(mode = "list", length = length(TIC_align)) # create empty list for binned data

names(TIC_binned) <- names(TIC_align) # copy names



# Perform binning for all TICs

for (q in 1:length(TIC_align)) {

  df <- TIC_align[[q]]

  RT_min <- round(min(TIC_align[[q]]$time_corrected), digits = 2) # determine minimum RT

  RT_max <- round(max(TIC_align[[q]]$time_corrected), digits = 2) # determine maximum RT

  breaks <- seq(from = RT_min, to = RT_max, by = delta) # define breaks between min and max

  time_bin <- seq(from = (RT_min + delta / 2), to = (RT_max), by = delta) # define bins

  intensity_bin <- matrixStats::binMeans(df$intensity, x = df$time_corrected, bx = breaks) # compute sample means in non-overlapping bins

  intensity_bin[is.nan(intensity_bin)] <- 0 # set missing values to 0

  binned <- data.frame(time_bin, intensity_bin) # create data frame with results

  TIC_binned[[q]] <- binned # add data frame to list

  rm(df, RT_min, RT_max, breaks, time_bin, intensity_bin, binned) # remove variables that we don't need anymore

} # close for loop

rm(delta) # remove variables that we don't need anymore



# Save object

saveRDS(object = TIC_binned, file = "TIC_validation_KRetI_binned.RDS")

Chapter 5 - Normalization

Normalization can be performed in different ways, two options are presented here: standard normal variate (SNV) and multiplicative scatter correction (MSC). Select the one that fits your data best.

SNV normalization performs normalization per individual TIC

TIC_normalized <- vector(mode = "list", length = length(TIC_binned)) # create empty list for normalized data

names(TIC_normalized) <- names(TIC_binned) # copy names



# Perform normalization using SNV

for (q in 1:length(TIC_binned)) {

  df <- TIC_binned[[q]]

  df$intensity_norm <- spectacles::snv(df$intensity_bin) # perform SNV normalization

  normalized <- data.frame(df$time_bin, df$intensity_norm)

  colnames(normalized) <- c("time", "intensity")

  TIC_normalized[[q]] <- normalized

  rm(df, normalized)

} # close for loop



# Change list into data matrix (required for further processing and PCA)

# Find minimum and maximum for each TIC

rt_min <- c()

rt_max <- c()

for (i in 1:length(TIC_normalized)) {

  temp_min <- min(TIC_normalized[[i]][["time"]])

  temp_max <- max(TIC_normalized[[i]][["time"]])

  rt_min <- c(rt_min, temp_min)

  rt_max <- c(rt_max, temp_max)

  rm(temp_min, temp_max)

} # close for loop

rm(i, q) # clean up environment



# Determine lowest maximum and highest minimum for the cutoff

rt_min <- max(rt_min)

rt_max <- min(rt_max)



# Define empty data frame

data.snv <- data.frame()



# Organise data into matrix (incl. RT cut-off at rt_min and rt_max)

for (j in 1:length(TIC_normalized)) {

  clipb <- TIC_binned[[j]]

  selection <- which(round(clipb$time, digits = 3) >= rt_min & round(clipb$time, digits = 3) <= rt_max)

  clipb <- clipb[selection, ]

  clipb <- as.data.frame(t(clipb))

  colnames(clipb) <- as.character(clipb["time", ])

  clipb <- clipb[-1, ]

  rownames(clipb) <- names(TIC_normalized[j])

  data.snv <- rbind(data.snv, clipb)

  rm(clipb)

} # close for loop

rm(j, rt_max, rt_min, TIC_normalized) # clean up environment



# Save object

saveRDS(object = data.snv, file = "TIC_validation_KRetI_binned_SNV.RDS")

MSC normalization performs normalization alongside a reference chromatogram

# Create data matrix which is required for MSC

# Find highest minimum and lowest maximum for the cutoff (this is done to be able to make a filled matrix for the PCA)

rt_min <- c()

rt_max <- c()



# Find minimum and maximum for each TIC

for (i in 1:length(TIC_binned)) {

  temp_min <- min(TIC_binned[[i]][["time_bin"]])

  temp_max <- max(TIC_binned[[i]][["time_bin"]])

  rt_min <- c(rt_min, temp_min)

  rt_max <- c(rt_max, temp_max)

  rm(temp_min, temp_max)

}

rm(i)



# Determine lowest maximum and highest minimum for the cutoff

rt_min <- max(rt_min)

rt_max <- min(rt_max)



# Define empty data frame

data.msc <- data.frame()



# Organise data for MSC (incl. RT cut-off at rt_min and rt_max)

for (j in 1:length(TIC_binned)) {

  clipb <- TIC_binned[[j]]

  selection <- which(round(clipb$time, digits = 3) >= rt_min & round(clipb$time, digits = 3) <= rt_max)

  clipb <- clipb[selection, ]

  clipb <- as.data.frame(t(clipb))

  colnames(clipb) <- as.character(clipb["time", ])

  clipb <- clipb[-1, ]

  rownames(clipb) <- names(TIC_binned[j])

  data.msc <- rbind(data.msc, clipb)

  rm(clipb)

} # close for loop

rm(j, rt_max, rt_min)

ref <- as.numeric(data.msc[ref.chromatogram, ]) # get values of reference chromatogram



# Perform normalization using MSC

TIC_normalized <- as.data.frame(pls::msc(X = as.matrix(data.msc), reference = ref))



# Save object

saveRDS(object = TIC_normalized, file = "TIC_validation_KRetI_binned_MSC.RDS")

Chapter 6 - Smoothing

Smoothing can be done in different ways as well. Three options are presented here: Savitzky-Golay smoothing, modified polynomial fitting and local windows & Gaussian weighting. Select the one that fits your data best.

Savitzky-Golay smoothing

# apply SGolay smoothing

TIC_smoothed <- as.data.frame(t(apply(TIC_normalized, 1, signal::sgolayfilt)))

colnames(TIC_smoothed) <- colnames(TIC_normalized) # add chromatogram names



# Save object

saveRDS(object = TIC_smoothed, file = "TIC_validation_KRetI_binned_SGolay.RDS")

Modified polynomial fitting

# Apply modified polynomial fitting

polyfit <- as.data.frame(baseline::baseline.modpolyfit(spectra = as.matrix(TIC_normalized))[[2]])

colnames(polyfit) <- colnames(TIC_normalized)

rownames(polyfit) <- rownames(TIC_normalized)



# Save object

saveRDS(file = "TIC_validation_KRetI_binned_polynomial.RDS", object = polyfit)

Local windows & Gaussian weighting

# Apply local window & Gaussian weighting

locsmooth <- as.data.frame(baseline::baseline.medianWindow(spectra = as.matrix(TIC_normalized), hwm = 0.05)[[2]])



# Save object

saveRDS(file = "TIC_validation_KRetI_binned_localwindow.RDS", object = locsmooth)

Chapter 7 - Derivatives

Another option is to calculate 1st and 2nd derivatives of the data. The code is presented here, you can test whether it is beneficial for your data.

First derivative

# Calculate first derivative

d1 <- as.data.frame(t(diff(t(TIC_smoothed), differences = 1)))



# Save object

saveRDS(object = as.data.frame(d1), file = "TIC_validation_KRetI_binned_1stderiv.RDS")

Second derivative

# Calculate second derivative

d2 <- as.data.frame(t(diff(t(TIC_smoothed), differences = 2)))



# Save object

saveRDS(object = as.data.frame(d2), file = "TIC_validation_KRetI_binned_2ndderiv.RDS")

Chapter 8 - Principal Component Analysis (PCA)

A PCA can be used to reduce the dimensions of the data. The code is presented here. First, the data is prepared by removing everything before the solvent peak (t0) and removing the internal standard peaks. If information from the blank samples is available, impurities must be removed as well.

Data preparation

# Load data

data.pca <- readRDS("TIC_validation_KRetI_binned_MSC_1month_phenol.RDS")



# Remove everything before t0.cutoff (to be defined by the operator) min to exclude t0

t0.cutoff <- 4.465

data.pca <- data.pca[, -which(as.numeric(colnames(data.pca)) < t0.cutoff)]



# Remove internal standard peaks (range of 0.2 min around APEX)

## Define reference chromatogram (same as for KRetI alignment)

ref <- getRt(all.ref.IS = all.ref.IS, chromatogram = ref.chromatogram)

range <- 0.2



# Remove IS peaks

data.pca <- data.pca[, -which(as.numeric(colnames(data.pca)) > (ref$`chlorobenzene-d5`[1] - range) &

  as.numeric(colnames(data.pca)) < (ref$`chlorobenzene-d5`[1] + range))] # Remove chlorobenzene-d5

data.pca <- data.pca[, -which(as.numeric(colnames(data.pca)) > (ref$`dichlorobenzene-d4`[1] - range) &

  as.numeric(colnames(data.pca)) < (ref$`dichlorobenzene-d4`[1] + range))] # Remove 1,4-dichlorobenzene-d4

data.pca <- data.pca[, -which(as.numeric(colnames(data.pca)) > (ref$`naphthalene-d8`[1] - range) &

  as.numeric(colnames(data.pca)) < (ref$`naphthalene-d8`[1] + range))] # Remove naphthalene-d8

data.pca <- data.pca[, -which(as.numeric(colnames(data.pca)) > (ref$`dibromobenzene-d4`[1] - range) &

  as.numeric(colnames(data.pca)) < (ref$`dibromobenzene-d4`[1] + range))] # Remove 1,4-dibromobenzene-d4

data.pca <- data.pca[, -which(as.numeric(colnames(data.pca)) > (ref$`phenanthrene-d10`[1] - range) &

  as.numeric(colnames(data.pca)) < (ref$`phenanthrene-d10`[1] + range))] # Remove phenanthrene-d10

data.pca <- data.pca[, -which(as.numeric(colnames(data.pca)) > (ref$`terbuthylazine-d5`[1] - range) &

  as.numeric(colnames(data.pca)) < (ref$`terbuthylazine-d5`[1] + range))] # Remove terbuthylazine-d5

Perform PCA

pc <- prcomp(data.pca, scale. = T) # perform PCA



data.pca$class <- sapply(rownames(data.pca), getClassification) # add class to data

data.pca$year <- sapply(rownames(data.pca), getYear) # add year to data

Create scree plot

factoextra::fviz_eig(pc, addlabels = TRUE)

Create biplot

factoextra::fviz_pca_biplot(pc,

  select.var = list(contrib = 15),

  col.var = "black", alpha.var = 0.2,

  repel = FALSE,

  habillage = data.pca$year,

  ggtheme = theme_minimal(),

  title = "PCA biplot",

  labelsize = 2

)

Chapter 9 - Hierarchical Clustering Analysis (HCA)

HCA can be used to generate clusters in the data. The code is presented here. First, the data is prepared by removing everything before the solvent peak (t0) and removing the internal standard peaks. If information from the blank samples is available, impurities must be removed as well. User input required in the section below, he/she can define the t0.cutoff and remove other irrelevant peaks if desired.

Data preparation

# Load data

data.hca <- readRDS("TIC_validation_KRetI_binned_MSC_1month_phenol.RDS")



# Remove everything before t0.cutoff (to be defined by the operator) min to exclude t0

t0.cutoff <- 4.465

data.hca <- data.hca[, -which(as.numeric(colnames(data.hca)) < t0.cutoff)]



# Remove internal standard peaks (range of 0.2 min around APEX)

## Define reference chromatogram (same as for KRetI alignment)

ref <- getRt(all.ref.IS = all.ref.IS, chromatogram = ref.chromatogram)

range <- 0.2



# Remove IS peaks

data.hca <- data.hca[, -which(as.numeric(colnames(data.hca)) > (ref$`chlorobenzene-d5`[1] - range) &

  as.numeric(colnames(data.hca)) < (ref$`chlorobenzene-d5`[1] + range))] # Remove chlorobenzene-d5

data.hca <- data.hca[, -which(as.numeric(colnames(data.hca)) > (ref$`dichlorobenzene-d4`[1] - range) &

  as.numeric(colnames(data.hca)) < (ref$`dichlorobenzene-d4`[1] + range))] # Remove 1,4-dichlorobenzene-d4

data.hca <- data.hca[, -which(as.numeric(colnames(data.hca)) > (ref$`naphthalene-d8`[1] - range) &

  as.numeric(colnames(data.hca)) < (ref$`naphthalene-d8`[1] + range))] # Remove naphthalene-d8

data.hca <- data.hca[, -which(as.numeric(colnames(data.hca)) > (ref$`dibromobenzene-d4`[1] - range) &

  as.numeric(colnames(data.hca)) < (ref$`dibromobenzene-d4`[1] + range))] # Remove 1,4-dibromobenzene-d4

data.hca <- data.hca[, -which(as.numeric(colnames(data.hca)) > (ref$`phenanthrene-d10`[1] - range) &

  as.numeric(colnames(data.hca)) < (ref$`phenanthrene-d10`[1] + range))] # Remove phenanthrene-d10

data.hca <- data.hca[, -which(as.numeric(colnames(data.hca)) > (ref$`terbuthylazine-d5`[1] - range) &

  as.numeric(colnames(data.hca)) < (ref$`terbuthylazine-d5`[1] + range))] # Remove terbuthylazine-d5



# According to Schollee et al. 2018: divide by maximum values

max_row_pos <- apply(data.hca, 1, max)

normmax.pos <- data.hca[, 1:ncol(data.hca)] / max_row_pos



# If desired, values with low intensities (below threshold) can be turned into 0, this can also be done the other way around (exclude high intensities)

# treshold <- 0.3

# normmax.pos[normmax.pos < threshold] <- 0

Perform HCA

pheatmap::pheatmap(normmax.pos,

  scale = "none",

  clustering_distance_rows = "euclidean",

  show_rownames = T,

  border_color = NA,

  cluster_cols = FALSE,

  main = "Hierarchical clustering of GC-MS analyses (Euclidean distance, max normalized)",

  fontsize_col = 4,

  width = 35,

  height = 20

)

The plot generated here is too small. It can be zoomed in when saved to a .png file, which can be done by simply adding a filename to the pheatmap() function.

Chapter 10 - Retrieve scan of interest

When a peak of interest has been found, the corresponding MS scan can be retrieved using the following code. Here, a .txt file will be generated which can be uploaded on MassBank Europe and MassBank of North America. User input required in the section below, he/she must select the chromatogram and scan of interest.

Get scan of interest

# Define chromatogram of interest

chrom.of.interest <- "180911_LOB_18_SS" # To be filled in by the operator



# Define (corrected) retention time of interest, retrieved from HCA plot

rt.of.interest <- 5.494 # To be filled in by the operator



# Find original retention time which belongs to the rt.of.interest

rt.original <- TIC_align[[chrom.of.interest]]$time[which.min(abs(TIC_align[[chrom.of.interest]]$time_corrected - rt.of.interest))]



# Convert retention time from minutes to seconds

rt.original.sec <- rt.original * 60



# Find MS spectrum which belongs to this retention time

# Retrieve all scans and their corresponding retention times in seconds. This will be used for the Rt window to look for the IS

all.Rt.scans <- NULL

for (y in 1:length(spectra[[chrom.of.interest]])) {

  Rt.scans <- data.frame(rtinseconds = spectra[[chrom.of.interest]][[y]]$rtinseconds, scan = spectra[[chrom.of.interest]][[y]]$scan)

  all.Rt.scans <- rbind(all.Rt.scans, Rt.scans) # append the scans and Rt to a vector

}



# Retrieve scan number of interest

scan.of.interest <- all.Rt.scans$scan[which.min(abs(all.Rt.scans$rtinseconds - rt.original.sec))]



# Retrieve mass spectrum

MS <- data.frame(spectra[[chrom.of.interest]][[scan.of.interest]][["mZ"]], spectra[[chrom.of.interest]][[scan.of.interest]][["intensity"]])



# Make intensities relative and clean up data frame

MS$relativeIntensity <- MS$spectra..chrom.of.interest....scan.of.interest.....intensity... / max(MS$spectra..chrom.of.interest....scan.of.interest.....intensity...) * 999

MS$spectra..chrom.of.interest....scan.of.interest.....intensity... <- NULL

colnames(MS) <- c("mZ", "relativeIntensity")

Export scan to .txt ready for uploading to MassBank Europe and MoNA

write.table(MS, file = paste0(chrom.of.interest, "_scan", scan.of.interest, ".txt"), row.names = FALSE)



# Export only top 20 peaks

MS.top20 <- data.table::setorder(MS, -relativeIntensity)[1:20, ]

write.table(MS.top20, file = paste0(chrom.of.interest, "_scan", scan.of.interest, "_top20.txt"), row.names = FALSE)

##          mZ relativeIntensity

## 1  94.03635         999.00000

## 2  66.08657         501.47166

## 3  65.08504         324.92512

## 4  39.07484         218.90578

## 5  40.08612         112.41123

## 6  55.06335          92.58730

## 7  63.05878          82.65805

## 8  95.07397          70.97227

## 9  38.06773          69.10941

## 10 43.06993          61.67833

## 11 50.06155          58.62468

## 12 51.06593          52.08628

## 13 62.06422          47.40763

## 14 41.09200          40.79415

## 15 67.08057          38.65261

## 16 53.05008          34.30169

## 17 42.08316          33.90247

## 18 64.08242          33.35000

## 19 47.07994          33.05748

## 20 61.06101          32.81697