403817 - LC-HRMS data analysis

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

When user input is necessary, it is represented like this: user input required. Various functions from patRoon 2.0 are applied in this script, more information on patRoon can be found in the tutorial, the handbook and in Helmus et al. 2021.

Initialization

First, required packages must be loaded.

# Load required packages

library(patRoon) # for data pre-processing and analysis - make sure patRoon 2.0 is installed!

library(GGally) # for plotting trends

library(ggplot2) # for plotting in general

library(viridis) # for color scales suitable for color-blind people

library(xlsx) # for writing results to an Excel file

Load data

Every folder contains samples from 3 locations and these have been measured in triplicate, the blank has been measured in triplicate as well. Sometimes exceptions occur (failed measurements), these are corrected.

First, the files in every folder are organised in a data frame per folder. This must be adjusted by the user for the dataset of interest.

anaInfoPos1 <- generateAnalysisInfo(

  paths = "\\\\nwg\\dfs\\projecttemp\\P403817_001\\HWL data calibrated\\25-3-2019\\pos",

  groups = c(rep("IJM-pos-01", 3), rep("MAAS-pos-01", 3), rep("LEKKAN-RIJN-pos-01", 3), rep("PROCBL-pos-01", 3)),

  blanks = rep("PROCBL-pos-01", 12)

)



anaInfoPos2 <- generateAnalysisInfo(

  paths = "\\\\nwg\\dfs\\projecttemp\\P403817_001\\HWL data calibrated\\23-4-2019\\pos",

  groups = c(rep("IJM-pos-02", 3), rep("MAAS-pos-02", 3), rep("LEKKAN-RIJN-pos-02", 3), rep("PROCBL-pos-02", 4)),

  blanks = rep("PROCBL-pos-02", 13)

)

# remove "PROCBL_GA7_01_6074.d" since this measurement contains no data

anaInfoPos2 <- anaInfoPos2[-which(anaInfoPos2$analysis == "PROCBL_GA7_01_6074"), ]

row.names(anaInfoPos2) <- NULL



# The files in this folder cause a fatal error during findFeatures(), so these measurements have been excluded for the analysis.

# anaInfoPos3 <- generateAnalysisInfo(

#   paths = "\\\\nwg\\dfs\\projecttemp\\P403817_001\\HWL data calibrated\\17-6-2019\\pos",

#   groups = c(rep("IJM-pos-03", 3), rep("MAAS-pos-03", 3), rep("LEKKAN-RIJN-pos-03", 3), rep("PROCBL-pos-03", 3)),

#   blanks = rep("PROCBL-pos-03", 12)

# )



anaInfoPos4 <- generateAnalysisInfo(

  paths = "\\\\nwg\\dfs\\projecttemp\\P403817_001\\HWL data calibrated\\15-7-2019\\pos",

  groups = c(rep("IJM-pos-04", 3), rep("MAAS-pos-04", 3), rep("LEKKAN-RIJN-pos-04", 3), rep("PROCBL-pos-04", 3)),

  blanks = rep("PROCBL-pos-04", 12)

)



anaInfoPos5 <- generateAnalysisInfo(

  paths = "\\\\nwg\\dfs\\projecttemp\\P403817_001\\HWL data calibrated\\12-8-2019\\pos",

  groups = c(rep("IJM-pos-05", 3), rep("MAAS-pos-05", 3), rep("LEKKAN-RIJN-pos-05", 3), rep("PROCBL-pos-05", 3)),

  blanks = rep("PROCBL-pos-05", 12)

)



anaInfoPos6 <- generateAnalysisInfo(

  paths = "\\\\nwg\\dfs\\projecttemp\\P403817_001\\HWL data calibrated\\9-9-2019\\pos",

  groups = c(rep("IJM-pos-06", 3), rep("MAAS-pos-06", 3), rep("LEKKAN-RIJN-pos-06", 3), rep("PROCBL-pos-06", 3)),

  blanks = rep("PROCBL-pos-06", 12)

)



anaInfoPos7 <- generateAnalysisInfo(

  paths = "\\\\nwg\\dfs\\projecttemp\\P403817_001\\HWL data calibrated\\7-10-2019\\pos",

  groups = c(rep("IJM-pos-07", 3), rep("MAAS-pos-07", 3), rep("LEKKAN-RIJN-pos-07", 3), rep("PROCBL-pos-07", 4)),

  blanks = rep("PROCBL-pos-07", 13)

)

# remove "PROCBL F_GA7_01_7630.d" since this measurement contains no data

anaInfoPos7 <- anaInfoPos7[-which(anaInfoPos7$analysis == "PROCBL F_GA7_01_7630"), ]

row.names(anaInfoPos7) <- NULL



anaInfoPos8 <- generateAnalysisInfo(

  paths = "\\\\nwg\\dfs\\projecttemp\\P403817_001\\HWL data calibrated\\4-11-2019\\pos",

  groups = c(rep("IJM-pos-08", 3), rep("MAAS-pos-08", 3), rep("LEKKAN-RIJN-pos-08", 3), rep("PROCBL-pos-08", 3)),

  blanks = rep("PROCBL-pos-08", 12)

)



anaInfoPos9 <- generateAnalysisInfo(

  paths = "\\\\nwg\\dfs\\projecttemp\\P403817_001\\HWL data calibrated\\2-12-2019\\pos",

  groups = c(rep("IJM-pos-09", 3), rep("MAAS-pos-09", 3), rep("LEKKAN-RIJN-pos-09", 3), rep("PROCBL-pos-09", 3)),

  blanks = rep("PROCBL-pos-09", 12)

)



anaInfoPos10 <- generateAnalysisInfo(

  paths = "\\\\nwg\\dfs\\projecttemp\\P403817_001\\HWL data calibrated\\6-1-2020\\pos",

  groups = c(rep("MAAS-pos-10", 3), rep("IJM-pos-10", 3), rep("LEKKAN-RIJN-pos-10", 3), rep("PROCBL-pos-10", 3)),

  blanks = rep("PROCBL-pos-10", 12)

)



anaInfoPos11 <- generateAnalysisInfo(

  paths = "\\\\nwg\\dfs\\projecttemp\\P403817_001\\HWL data calibrated\\3-2-2020\\pos",

  groups = c(rep("IJM-pos-11", 3), rep("MAAS-pos-11", 3), rep("LEKKAN-RIJN-pos-11", 3), rep("PROCBL-pos-11", 3)),

  blanks = rep("PROCBL-pos-11", 12)

)



anaInfoPos12 <- generateAnalysisInfo(

  paths = "\\\\nwg\\dfs\\projecttemp\\P403817_001\\HWL data calibrated\\2-3-2020\\pos",

  groups = c(rep("IJM-pos-12", 3), rep("MAAS-pos-12", 3), rep("LEKKAN-RIJN-pos-12", 3), rep("PROCBL-pos-12", 3)),

  blanks = rep("PROCBL-pos-12", 12)

)



anaInfoPos13 <- generateAnalysisInfo(

  paths = "\\\\nwg\\dfs\\projecttemp\\P403817_001\\HWL data calibrated\\30-3-2020\\pos",

  groups = c(rep("IJM-pos-13", 3), rep("MAAS-pos-13", 3), rep("LEKKAN-RIJN-pos-13", 3), rep("PROCBL-pos-13", 3)),

  blanks = rep("PROCBL-pos-13", 12)

)



anaInfoPos14 <- generateAnalysisInfo(

  paths = "\\\\nwg\\dfs\\projecttemp\\P403817_001\\HWL data calibrated\\28-4-2020\\pos",

  groups = c(rep("IJM-pos-14", 3), rep("MAAS-pos-14", 3), rep("LEKKAN-RIJN-pos-14", 3), rep("PROCBL-pos-14", 3)),

  blanks = rep("PROCBL-pos-14", 12)

)



anaInfoPos15 <- generateAnalysisInfo(

  paths = "\\\\nwg\\dfs\\projecttemp\\P403817_001\\HWL data calibrated\\25-5-2020\\pos",

  groups = c(rep("IJM-pos-15", 3), rep("MAAS-pos-15", 3), rep("LEKKAN-RIJN-pos-15", 3), rep("PROCBL-pos-15", 3)),

  blanks = rep("PROCBL-pos-15", 12)

)



anaInfoPos16 <- generateAnalysisInfo(

  paths = "\\\\nwg\\dfs\\projecttemp\\P403817_001\\HWL data calibrated\\22-6-2020\\pos",

  groups = c(rep("IJM-pos-16", 3), rep("MAAS-pos-16", 3), rep("LEKKAN-RIJN-pos-16", 3), rep("PROCBL-pos-16", 3)),

  blanks = rep("PROCBL-pos-16", 12)

)



anaInfoPos17 <- generateAnalysisInfo(

  paths = "\\\\nwg\\dfs\\projecttemp\\P403817_001\\HWL data calibrated\\20-7-2020\\pos",

  groups = c(rep("IJM-pos-17", 3), rep("MAAS-pos-17", 3), rep("LEKKAN-RIJN-pos-17", 3), rep("PROCBL-pos-17", 3)),

  blanks = rep("PROCBL-pos-17", 12)

)



anaInfoPos18 <- generateAnalysisInfo(

  paths = "\\\\nwg\\dfs\\projecttemp\\P403817_001\\HWL data calibrated\\17-8-2020\\pos",

  groups = c(rep("IJM-pos-18", 3), rep("MAAS-pos-18", 3), rep("LEKKAN-RIJN-pos-18", 3), rep("PROCBL-pos-18", 3)),

  blanks = rep("PROCBL-pos-18", 12)

)



anaInfoPos19 <- generateAnalysisInfo(

  paths = "\\\\nwg\\dfs\\projecttemp\\P403817_001\\HWL data calibrated\\14-9-2020\\pos",

  groups = c(rep("IJM-pos-19", 3), rep("MAAS-pos-19", 3), rep("LEKKAN-RIJN-pos-19", 3), rep("PROCBL-pos-19", 3)),

  blanks = rep("PROCBL-pos-19", 12)

)



anaInfoPos20 <- generateAnalysisInfo(

  paths = "\\\\nwg\\dfs\\projecttemp\\P403817_001\\HWL data calibrated\\12-10-2020\\pos",

  groups = c(rep("IJM-pos-20", 3), rep("MAAS-pos-20", 3), rep("LEKKAN-RIJN-pos-20", 3), rep("PROCBL-pos-20", 3)),

  blanks = rep("PROCBL-pos-20", 12)

)



anaInfoPos21 <- generateAnalysisInfo(

  paths = "\\\\nwg\\dfs\\projecttemp\\P403817_001\\HWL data calibrated\\9-11-2020\\pos",

  groups = c(rep("IJM-pos-21", 3), rep("MAAS-pos-21", 3), rep("LEKKAN-RIJN-pos-21", 3), rep("PROCBL-pos-21", 3)),

  blanks = rep("PROCBL-pos-21", 12)

)



anaInfoPos22 <- generateAnalysisInfo(

  paths = "\\\\nwg\\dfs\\projecttemp\\P403817_001\\HWL data calibrated\\7-12-2020\\pos",

  groups = c(rep("IJM-pos-22", 3), rep("MAAS-pos-22", 3), rep("LEKKAN-RIJN-pos-22", 3), rep("PROCBL-pos-22", 3)),

  blanks = rep("PROCBL-pos-22", 12)

)



# Generate 1 data frame with all analyses and clean up environment

anaInfoPos <- rbind(

  anaInfoPos1, anaInfoPos2, anaInfoPos4, anaInfoPos5, anaInfoPos6, anaInfoPos7, anaInfoPos8, anaInfoPos9, anaInfoPos10, anaInfoPos11, anaInfoPos12, anaInfoPos13, anaInfoPos14, anaInfoPos15, anaInfoPos16, anaInfoPos17, anaInfoPos18, anaInfoPos19, anaInfoPos20, anaInfoPos21, anaInfoPos22

)



# clean up environment

rm(

  anaInfoPos1, anaInfoPos2, anaInfoPos4, anaInfoPos5, anaInfoPos6, anaInfoPos7, anaInfoPos8, anaInfoPos9, anaInfoPos10, anaInfoPos11, anaInfoPos12, anaInfoPos13, anaInfoPos14, anaInfoPos15, anaInfoPos16, anaInfoPos17, anaInfoPos18, anaInfoPos19, anaInfoPos20, anaInfoPos21, anaInfoPos22

)

row.names(anaInfoPos) <- NULL

This results in a data frame with the directory path of he file containing the analysis data, the name of the analysis (which is the file name without file extension), to which replicate group the analysis belongs and which replicate group should be used for blank subtraction.

Convert files

Next, the files are converted from the vendor format (Bruker) into an open source format (mzML) using the ProteoWizard algorithm.

doDataPretreatment <- TRUE # set to FALSE to exclude data pre-treatment

if (doDataPretreatment) {

  convertMSFiles(anaInfo = anaInfoPos, from = "bruker", to = "mzML", algorithm = "pwiz", centroid = TRUE)

}

Find features

All files are screened for features using the OpenMS algorithm and its default settings. This may take some processing time, especially for larger datasets. These parameters can be optimized by the user.

fList <- findFeatures(anaInfoPos, "openms",

  noiseThrInt = 1000, chromSNR = 3, chromFWHM = 5, mzPPM = 10,

  reEstimateMTSD = TRUE, traceTermCriterion = "sample_rate", traceTermOutliers = 5, minSampleRate = 0.5, minTraceLength = 3,

  maxTraceLength = -1, widthFiltering = "fixed", minFWHM = 1, maxFWHM = 30, traceSNRFiltering = FALSE, localRTRange = 10,

  localMZRange = 6.5, isotopeFilteringModel = "metabolites (5% RMS)", MZScoring13C = FALSE, useSmoothedInts = TRUE, extraOpts = NULL,

  intSearchRTWindow = 3, useFFMIntensities = FALSE, verbose = TRUE

)

All features before t0 (the solvent peak) are removed prior to grouping, this can be adjusted by the user.

lower.rt <- 180 # everything before this Rt (in seconds) will be removed

upper.rt <- Inf # everything after this Rt (in seconds) will be removed

fList <- filter(fList, retentionRange = c(lower.rt, upper.rt))

## Done! Filtered 18094 (5.38%) features. Remaining: 318009

The features are grouped and aligned between analyses. Used OpenMS algorithm here with all default options except maxGroupRT, increased to 30, so that internal standard atrazin-d5 gets grouped correctly. These parameters can be adjusted by the user.

fGroups <- groupFeatures(fList, "openms",

  rtalign = TRUE, QT = FALSE, maxAlignRT = 30,

  maxAlignMZ = 0.005, maxGroupRT = 30, maxGroupMZ = 0.005, extraOptsRT = NULL,

  extraOptsGroup = NULL, verbose = TRUE

)

Apply rule-based filtering: only remove features that aren’t in all analysis of replicate groups and have intensity RSD (relative standard deviation) >75%. These parameters can be adjusted by the user.

fGroups <- filter(fGroups,

  preAbsMinIntensity = 100, absMinIntensity = 0, relMinReplicateAbundance = 1,

  maxReplicateIntRSD = 0.75, blankThreshold = 0, removeBlanks = FALSE,

  retentionRange = NULL, mzRange = NULL

)

## Applying intensity filter... Done! Filtered 0 (0.00%) features and 0 (0.00%) feature groups. Remaining: 242541 features in 9398 groups.

## Applying replicate abundance filter... Done! Filtered 0 (0.00%) features and 0 (0.00%) feature groups. Remaining: 242541 features in 9398 groups.

## Applying blank filter... Done! Filtered 0 (0.00%) features and 0 (0.00%) feature groups. Remaining: 242541 features in 9398 groups.

Normalization

The feature intensities must be normalized according to the intensity of the internal standard. Normalization is performed by dividing feature intensities by the intensity of the internal standard, in this case, atrazine-d5.

The user needs to define the feature group that corresponds to the internal standard here.

# find and enter feature group that represents atrazin-d5

rtmin <- 540 # lower Rt limit where atrazin-d5 is expected, in seconds

rtmax <- 570 # upper Rt limit where atrazin-d5 is expected, in seconds

mzmin <- 220 # lower m/z limit of atrazin-d5

mzmax <- 222 # upper m/z limit of atrazin-d5

temp <- filter(fGroups, retentionRange = c(rtmin, rtmax), mzRange = c(mzmin, mzmax))

## Applying retention filter... Done! Filtered 236010 (97.31%) features and 9063 (96.44%) feature groups. Remaining: 6531 features in 335 groups.

## Applying mz filter... Done! Filtered 6270 (96.00%) features and 332 (99.10%) feature groups. Remaining: 261 features in 3 groups.

# Choose correct feature group from 'temp'

fgroup.is <- "M221_R550_5692"

# create data frame

fGroups.df <- as.data.frame(fGroups)



# normalize data using the intensity of internal standard atrazin-d5

ind <- which(fGroups.df$group == fgroup.is) # find atrazin-d5 in total set

is <- as.numeric(fGroups.df[ind, -c(1:3)]) # get intensities of atrazin-d5 (and remove columns 1:3 containing feature information (mz, Rt))

normalized <- as.data.frame(mapply("/", fGroups.df[, -c(1:3)], is)) # divide every intensity by atrazin-d5 intensity of the same measurement

normalized <- normalized * 10000 # multiply all values with 10,000 in order to get whole numbers

normalized$fGroups <- fGroups.df$group # add feature group column again since it was lost during mapply

fGroups.df <- normalized

rm(ind, is, normalized, rtmin, rtmax, mzmin, mzmax, temp) # clean up environment

Intensity filter

An intensity threshold can be set to remove all features with intensities below this value, this threshold is relative to the internal standard intensity. The user can adjust this threshold.

# due to the normalization step, the threshold is the same everywhere

threshold.pct <- 10 # define intensity percentage for the threshold

threshold <- (threshold.pct / 100) * 10000 # calculate absolute value



# apply threshold, set everything to zero if it is below the threshold

fGroups.df[fGroups.df < threshold.pct] <- 0



# clean up environment

rm(threshold.pct, threshold)

Re-order data

Create per sample a column with average feature intensities.

# get unique sample groups

groups <- unique(anaInfoPos$group)



# get unique sample group & blank combinations

groups.blanks <- unique(anaInfoPos[, c("group", "blank")])

row.names(groups.blanks) <- NULL



# create empty data frame and add feature groups as an extra column

fGroups.averaged <- setNames(data.frame(matrix(ncol = length(groups), nrow = nrow(fGroups.df))), groups)

fGroups.averaged$feature <- fGroups.df$fGroups # add column with feature groups



# fill data frame for every sample with the average intensity of the three replicates

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

  temp <- anaInfoPos$analysis[which(anaInfoPos$group == groups[i])]

  fGroups.averaged[, groups[i]] <- rowMeans(fGroups.df[, temp])

  rm(temp)

}

rm(i)

head(fGroups.averaged[, c(1:5, ncol(fGroups.averaged))])

##   IJM-pos-01 MAAS-pos-01 LEKKAN-RIJN-pos-01 PROCBL-pos-01 IJM-pos-02

## 1      0.000        0.00                  0             0      0.000

## 2      0.000        0.00                  0             0      0.000

## 3   2869.443     3137.47                  0             0   3936.567

## 4      0.000        0.00                  0             0      0.000

## 5      0.000        0.00                  0             0      0.000

## 6      0.000        0.00                  0             0      0.000

##       feature

## 1 M39_R1014_1

## 2  M39_R923_2

## 3  M39_R954_3

## 4  M39_R884_4

## 5  M43_R283_8

## 6 M43_R294_10

In some cases, the internal standard cannot be detected. These samples must be removed prior to perform blank correction. The user needs to modify the removal of measurements.

# create new data frame for the blank corrected data

fGroups.corrected <- fGroups.averaged



# check if there are any columns filled with NaN (due to the internal standard that could not be found)

toremove <- c()

for (i in 1:length(colnames(fGroups.corrected))) {

  x <- !any(is.nan(fGroups.corrected[, i]))

  if (x == FALSE) {

    toremove <- c(toremove, i) # remove column if NaN's are present

  }

}



# inspect filetype that has to be removed

colnames(fGroups.corrected)[toremove]

## [1] "PROCBL-pos-07"

rm(toremove, x, i) # clean up environment



# remove all measurements with pos-07 - This is an example, can be adjusted by the user

fGroups.corrected <- fGroups.corrected[, -(base::grep(pattern = "-pos-07", x = colnames(fGroups.corrected)))]

groups <- groups[-(base::grep(pattern = "-pos-07", x = groups))]

groups.blanks <- groups.blanks[-(base::grep(pattern = "-pos-07", x = groups.blanks$group)), ]

row.names(groups.blanks) <- NULL

Blank correction

All feature groups with intensities < 5x intensity in the blank samples are removed. Similarly to the previous step, this is a common step that allows to reduce the number of features to be analysed by removing features which are present in blanks in substantial levels. Finally, 1x the intensity of the blank is removed from features which are also present in blanks, although with an intensity of at least 5 times compared to the blank. This is done in order to correct for the variation in background between batches and measurement series.

# check if feature exceeds 5x blank intensity. if so, subtract blank intensity.

# if feature does not exceed 5x blank intensity, turn intensity into 0.

for (k in 1:length(groups)) {

  blank <- groups.blanks$blank[which(groups.blanks$group == groups[k])] # get correct blank column

  temp <- fGroups.averaged[, groups[k]] - (5 * fGroups.averaged[, blank]) # subract 5x blank intensity

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

    if (temp[n] <= 0) { # if feature intensity is below or equal to 0, it does not exceed 5x blank intensity so turn intensity into (or remain) 0

      fGroups.corrected[n, groups[k]] <- 0

    } else { # if feature intensity is > 0, it exceeds 5x blank intensity, so it should remain but 1x blank intensity must be subtracted.

      fGroups.corrected[n, groups[k]] <- fGroups.averaged[n, groups[k]] - fGroups.averaged[n, blank]

    }

  }

  rm(blank, temp, n)

}



# remove columns representing the blank measurements

fGroups.corrected <- fGroups.corrected[, -(base::grep(pattern = "PROCBL", x = colnames(fGroups.corrected)))]



# extra clean up set to remove columns with NaN (occurs when internal standard cannot be found)

tokeep <- c()

for (i in 1:length(colnames(fGroups.corrected))) {

  x <- !any(is.nan(fGroups.corrected[, i]))

  if (x == TRUE) {

    tokeep <- c(tokeep, i) # only keep column if no NaN's are present

  }

}

fGroups.corrected <- fGroups.corrected[, tokeep]



# remove feature groups that are not present in any of the measurements

ind.featureID <- which(colnames(fGroups.corrected) == "feature")

fGroups.corrected <- fGroups.corrected[apply(fGroups.corrected[, -ind.featureID], 1, function(x) !all(x == 0)), ]

row.names(fGroups.corrected) <- NULL



# clean up environment

rm(groups, i, ind.featureID, k, x, groups.blanks, fGroups.df, fGroups.averaged, tokeep)



head(fGroups.corrected[, c(1:5, ncol(fGroups.corrected))])

##   IJM-pos-01 MAAS-pos-01 LEKKAN-RIJN-pos-01 IJM-pos-02 MAAS-pos-02     feature

## 1      0.000        0.00                  0      0.000           0 M39_R1014_1

## 2      0.000        0.00                  0      0.000           0  M39_R923_2

## 3   2869.443     3137.47                  0   3936.567           0  M39_R954_3

## 4      0.000        0.00                  0      0.000           0 M43_R294_10

## 5      0.000        0.00                  0      0.000           0 M43_R288_13

## 6      0.000        0.00                  0      0.000           0 M43_R424_14

Subset location

In this case a subset of the samples from the river Meuse is taken. The user can adjust this section to select a subset of interest.

# Select only data from the river Meuse

x.subset <- base::grep(pattern = "MAAS", x = colnames(fGroups.corrected))

x.feature <- base::grep(pattern = "feature", x = colnames(fGroups.corrected))

fGroups.df.subset <- fGroups.corrected[, c(x.feature, x.subset)]



# remove features that are not present in any of the measurements selected for this subset

x.feature <- base::grep(pattern = "feature", x = colnames(fGroups.df.subset))

fGroups.df.subset <- fGroups.df.subset[apply(fGroups.df.subset[, -x.feature], 1, function(x) !all(x == 0)), ]

row.names(fGroups.df.subset) <- NULL



rm(x.subset, x.feature)

Hierarchical Cluster Analysis

Hierarchical Cluster Analyses can be used for anomaly detection, clusters and samples of interest can be selected easily. Here Pearson’s correlation and Ward.D2’s method for clustering are used. Because of the large variability in feature intensities, these are log-transformed prior to HCA.

# log10 transformation 

data.hca <- log10(fGroups.df.subset[, 2:ncol(fGroups.df.subset)] + 1)

row.names(data.hca) <- fGroups.df.subset$feature



# remove features with only one occurrence

toremove <- c()

for (i in 1:nrow(data.hca)) {

  intensities <- as.numeric(data.hca[i, ])

  zeros <- intensities == 0

  if (ncol(data.hca) - (sum(zeros, na.rm = TRUE)) <= 1) {

    toremove <- c(toremove, i)

  }

  rm(zeros, intensities)

}

data.hca <- data.hca[-toremove, ]

rm(i, toremove)



# add feature name

annotation.row.subset <- row.names(data.hca)



# clustering of features based on pearson correlation

result <- pheatmap::pheatmap(data.hca,

  scale = "none", show_rownames = F,

  labels_col = colnames(data.hca),

  cutree_rows = 9,

  cluster_cols = FALSE,

  clustering_distance_rows = "correlation", # Pearson correlation

  clustering_method = "ward.D2",

  main = "log10 transformed, ward D2, Pearson correlation",

  fontsize_row = 4,

  width = 35,

  height = 20

)

# add cluster annotation and information on retention time and m/z

temp.result <- data.frame(cutree(result$tree_row, k = 9))

colnames(temp.result) <- "Cluster"

temp.result$Cluster <- as.character(temp.result$Cluster)

temp.result$featureID <- rownames(temp.result)

info <- as.data.frame(fGroups)[, 1:3]

names(info)[names(info) == "group"] <- "featureID"

info$ret <- info$ret / 60

annotation.row.result2 <- merge(x = temp.result, y = info, by = "featureID", all.y = FALSE)

row.names(annotation.row.result2) <- annotation.row.result2$featureID

annotation.row.result2$featureID <- NULL

names(annotation.row.result2)[names(annotation.row.result2) == "ret"] <- "Rt"

rm(temp.result, info) # clean environment



# Define color scale

dat.breaks <- quantile(data.hca$`MAAS-pos-01`, probs = seq(0, 1, length.out = 40))

dat.breaks <- dat.breaks[!duplicated(dat.breaks)]



# Define color palette for mass, retention time and clusters 

ann.colors <- list(

  Rt = c("#004c6d", "#346888", "#5886a5", "#7aa6c2", "#9dc6e0", "#c1e7ff"),

  mz = c("#8e0152", "#a63d6e", "#bd648c", "#d48aaa", "#eaafc8", "#ffd5e7"),

  Cluster = c("1" = "#001219", "2" = "#005F73", "3" = "#0A9396", "4" = "#94D2BD", "5" = "#E9D8A6", "6" = "#EE9B00", "7" = "#CA6702", "8" = "#BB3E03", "9" = "#9B2226")

)



# clustering of features based on Pearson correlation including cluster categories

result <- pheatmap::pheatmap(data.hca,

  scale = "none", show_rownames = F,

  color = viridis::viridis(length(dat.breaks)),

  breaks = dat.breaks,

  labels_col = colnames(data.hca),

  cutree_rows = 9,

  cluster_cols = FALSE,

  clustering_distance_rows = "correlation", # Pearson correlation

  clustering_method = "ward.D2", # Ward's minimum variance method

  main = "log10 transformed, ward.D2, Pearson correlation",

  annotation_row = annotation.row.result2,

  annotation_colors = ann.colors,

  annotation_names_row = TRUE,

  annotation_legend = TRUE,

  fontsize_row = 4,

  width = 35,

  height = 20,

)

# Clean up environment

rm(ann.colors, annotation.row.result2, dat.breaks)



# Create data frame with cluster number per feature group

fGroups.clusters <- data.frame(cutree(result$tree_row, k = 9))

fGroups.clusters$feature <- row.names(fGroups.clusters)

row.names(fGroups.clusters) <- NULL

colnames(fGroups.clusters) <- c("Cluster", "FeatureID")

head(fGroups.clusters)

##   Cluster   FeatureID

## 1       1  M39_R954_3

## 2       2 M43_R288_13

## 3       1 M45_R222_28

## 4       3 M45_R270_29

## 5       2 M45_R919_35

## 6       1 M45_R961_39

Principal Component Analysis

Principal Component Analysis can be used to highlight the presence of groups within samples or to detect the presence of outliers. Because of the large number of features still present in the dataset and their variability, the amount of variance explained by the first principal components is very low in this example.

# Perform PCA on all data

data.pca <- as.data.frame(t(fGroups.corrected))

names(data.pca) <- data.pca["feature", ]

data.pca <- data.pca[-which(row.names(data.pca) == "feature"), ]

data.pca <- dplyr::mutate_all(data.pca, as.numeric)

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



# Create scree plot

sp <- factoextra::fviz_eig(pc, addlabels = TRUE)



# Add extra column for habillage with sample name

data.pca$class <- substr(row.names(data.pca), 1, nchar(row.names(data.pca)) - 7)



# Biplot of samples and features

factoextra::fviz_pca_biplot(pc,

  select.var = list(contrib = 15),

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

  repel = FALSE,

  habillage = data.pca$class,

  ggtheme = theme_minimal(),

  title = "PCA biplot",

  labelsize = 2

)

Trend analysis

First, define functions to apply linear regression, Spearman rank correlation coefficient and Mann-Kendall rank correlation coefficient.

# define function to get logarithmic model

getLogLm <- function(intensities) {

  x <- 1:length(intensities)

  y <- as.numeric(intensities)

  # only perform function if 3 or more data points are available

  zeros <- y > 0

  if (sum(zeros, na.rm = TRUE) > 2) {

    model <- lm(y ~ log(x)) # generate model

    pvalue <- round(summary(model)$coefficients[2, 4], digits = 3) # extract p-value of the slope

    slope <- round(summary(model)$coefficients[2, 1], digits = 3) # extract slope

    output <- paste(pvalue, slope, sep = ";") # create output

  } else { # if less than 3 data points available, return NA

    output <- NA

  }

  return(output)

}



# define function to get linear model

getLinLm <- function(intensities) {

  x <- 1:length(intensities)

  y <- as.numeric(intensities)

  # only perform function if 3 or more data points are available

  zeros <- y > 0

  if (sum(zeros, na.rm = TRUE) > 2) {

    model <- lm(y ~ x)

    pvalue <- round(summary(model)$coefficients[2, 4], digits = 3) # extract p-value of the slope

    slope <- round(summary(model)$coefficients[2, 1], digits = 3) # extract slope

    output <- paste(pvalue, slope, sep = ";") # create output

  } else { # if less than 3 data points available, return NA

    output <- NA

  }

  return(output)

}



# Spearman rank correlation coefficient

getSpearman <- function(intensities) {

  x <- 1:length(intensities)

  y <- as.numeric(intensities)

  # only perform function if 3 or more data points are available

  zeros <- y > 0

  if (sum(zeros, na.rm = TRUE) > 2) {

    temp <- cor.test(x, y, method = "spearman", exact = FALSE)

    pvalue <- round(temp$p.value, digits = 3)

    rho <- round(temp$estimate, digits = 3)

    output <- paste(pvalue, rho, sep = ";") # create output

  } else { # if less than 3 data points available, return NA

    output <- NA

  }

  return(output)

}



# Mann-Kendall correlation coefficient

getMannKendall <- function(intensities) {

  x <- 1:length(intensities)

  y <- as.numeric(intensities)

  # only perform function if 3 or more data points are available

  zeros <- y > 0

  if (sum(zeros, na.rm = TRUE) > 2) {

    temp <- cor.test(x, y, method = "kendall", exact = FALSE)

    pvalue <- round(temp$p.value, digits = 3)

    rho <- round(temp$estimate, digits = 3)

    output <- paste(pvalue, rho, sep = ";") # create output

  } else { # if less than 3 data points available, return NA

    output <- NA

  }

  return(output)

}

Next, apply the functions to the data set. The user can modify this to select a statistical test of interest.

# apply tests/models to the data, make sure that the correct intensity columns are selected

fGroups.df.subset$Loglm <- apply(fGroups.df.subset[, 2:21], 1, getLogLm)

fGroups.df.subset$Linlm <- apply(fGroups.df.subset[, 2:21], 1, getLinLm)

fGroups.df.subset$Spearman <- apply(fGroups.df.subset[, 2:21], 1, getSpearman)

fGroups.df.subset$MannKendall <- apply(fGroups.df.subset[, 2:21], 1, getMannKendall)



# split columns and transfer to a dataframe with a more general name

fGroups.df.trend <- tidyr::separate(data = fGroups.df.subset, col = Loglm, into = c("Loglm_pvalue", "Loglm_slope"), sep = ";")

fGroups.df.trend <- tidyr::separate(data = fGroups.df.trend, col = Linlm, into = c("Linlm_pvalue", "Linlm_slope"), sep = ";")

fGroups.df.trend <- tidyr::separate(data = fGroups.df.trend, col = Spearman, into = c("Spearman_pvalue", "Spearman_rho"), sep = ";")

fGroups.df.trend <- tidyr::separate(data = fGroups.df.trend, col = MannKendall, into = c("MannKendall_pvalue", "MannKendall_rho"), sep = ";")



# turn characters into numeric

fGroups.df.trend$Loglm_pvalue <- as.numeric(fGroups.df.trend$Loglm_pvalue)

fGroups.df.trend$Loglm_slope <- as.numeric(fGroups.df.trend$Loglm_slope)

fGroups.df.trend$Linlm_pvalue <- as.numeric(fGroups.df.trend$Linlm_pvalue)

fGroups.df.trend$Linlm_slope <- as.numeric(fGroups.df.trend$Linlm_slope)

fGroups.df.trend$Spearman_pvalue <- as.numeric(fGroups.df.trend$Spearman_pvalue)

fGroups.df.trend$Spearman_rho <- as.numeric(fGroups.df.trend$Spearman_rho)

fGroups.df.trend$MannKendall_pvalue <- as.numeric(fGroups.df.trend$MannKendall_pvalue)

fGroups.df.trend$MannKendall_rho <- as.numeric(fGroups.df.trend$MannKendall_rho)

Next, select the features that have a p-value < 0.05 and a positive rho value for the Mann-Kendall test.

# select features that show in Mann-Kendall a pvalue < 0.05

features.oi <- c()

for (i in 1:nrow(fGroups.df.trend)) {

  pvalue <- as.numeric(fGroups.df.trend[i, "MannKendall_pvalue"])

  if (is.na(pvalue)) {

    pvalue <- 50

  } else if (pvalue < 0.05) {

    features.oi <- c(features.oi, fGroups.df.trend[i, "feature"])

  }

  rm(pvalue)

}

features.oi <- fGroups.df.trend[which(fGroups.df.trend$feature %in% features.oi), ]

row.names(features.oi) <- NULL



# select features with a positive rho, indicating an increasing trend

features.incr <- c()

for (i in 1:nrow(features.oi)) {

  rho <- as.numeric(features.oi[i, "MannKendall_rho"])

  if (is.na(rho)) {

    rho <- 0

  } else if (rho > 0) {

    features.incr <- c(features.incr, features.oi[i, "feature"])

  }

}

features.incr <- features.oi[which(features.oi$feature %in% features.incr), ]

row.names(features.incr) <- NULL



# plot those features

plot <- GGally::ggparcoord(

  data = features.incr, columns = 2:21, groupColumn = 1, showPoints = TRUE,

  title = "Increasing feature groups", missing = "exclude", scale = "globalminmax"

) +

  theme_bw() +

  xlab("sample") +

  ylab("normalized intensity (AU)") +

  theme(

    plot.title = element_text(size = 10), axis.title.x = element_text(size = 8), axis.title.y = element_text(size = 8),

    legend.position = "none", axis.text.x = element_text(angle = 90, vjust = 0.5), legend.title = element_text(size = 8)

  )

plot

Annotation

Here, we try to tentatively identify a feature with an increasing trend. The user can define a feature of interest and sample location here.

foi <- "M120_R355_1744" # fill in feature of interest

location <- "MAAS" # fill in sample location

outfile <- paste0(foi, "_", location) # define filename

ind <- grep(pattern = location, anaInfoPos$group) # select all indices of the files of the location

fGroup.oi <- fGroups[ind, foi] # select feature group of interest, only from the location of interest



# Retrieve MS peak lists

avgMSListParams <- getDefAvgPListParams(clusterMzWindow = 0.002)

mslists <- generateMSPeakLists(fGroup.oi, "mzr",

  maxMSRtWindow = 5, precursorMzWindow = 4,

  avgFeatParams = avgMSListParams,

  avgFGroupParams = avgMSListParams

)

saveRDS(mslists, file = paste0(outfile, "_mslist.RDS")) # save mslists to .rds as it takes quite some computation time



# Rule based filtering of MS peak lists. You may want to tweak this. See the patRoon manual for more information.

mslists <- filter(mslists,

  absMSIntThr = NULL, absMSMSIntThr = NULL, relMSIntThr = NULL, relMSMSIntThr = 0.02,

  topMSPeaks = NULL, topMSMSPeaks = 10

)



# Calculate formula candidates

formulas.sirius <- generateFormulas(fGroup.oi, mslists, "sirius",

  relMzDev = 5, elements = "CHNOPSClFBr", profile = "qtof",

  calculateFeatures = TRUE, featThresholdAnn = 0.75, adduct = "[M+H]+"

)

formulas.genform <- generateFormulas(fGroup.oi, mslists, "genform",

  relMzDev = 5, elements = "CHNOPSClFBr",

  calculateFeatures = TRUE, featThresholdAnn = 0.75, adduct = "[M+H]+"

)



# Write results away to Excel file

xlsx::write.xlsx(as.data.frame(formulas.sirius), file = paste0(outfile, ".xlsx"), sheetName = "SIRIUS", row.names = FALSE)

xlsx::write.xlsx(as.data.frame(formulas.genform), file = paste0(outfile, ".xlsx"), sheetName = "GenForm", append = TRUE, row.names = FALSE)



# generate compounds for feature of interest

compounds <- generateCompounds(fGroup.oi, mslists, "metfrag",

  method = "CL", dbRelMzDev = 5, fragRelMzDev = 5, fragAbsMzDev = 0.002,

  adduct = "[M+H]+", database = "pubchemlite", maxCandidatesToStop = 2500

)



# add formula scoring to improve ranking

compounds <- addFormulaScoring(compounds, formulas.genform, updateScore = TRUE)



# add results to existing excel sheet

xlsx::write.xlsx(as.data.frame(compounds), file = paste0(outfile, ".xlsx"), sheetName = "MetFrag", append = TRUE, row.names = FALSE)

The MetFrag data can be used to assign fragments to the MS2 peaks. And the MS2 of the feature is compared here to an MS2 of benzotriazole from MassBankEU.

## Identifying 1 feature groups with MetFrag...

## Converting to algorithm specific adducts... Done!

## Loaded 26 compounds from 1 features (100.00%).