---
title: "Clustering Mass Spectra from Low Resolution GC-EI-MS Data Using CluMSID"
author: "Tobias Depke"
date: '`r format(Sys.Date(), "%B %d, %Y")`'
output: rmarkdown::html_vignette
vignette: >
    %\VignetteIndexEntry{CluMSID GC-EI-MS Tutorial}
    %\VignetteEngine{knitr::rmarkdown}
    %\VignetteEncoding{UTF-8}
    %\VignetteDepends{metaMS, metaMSdata, CluMSIDdata, xcms}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>",
    echo = TRUE
)
```

<style>
p.caption {
    font-size: 80%;
    text-align: left;
}
</style>

```{r captions, include=FALSE}
fig1 <- paste("**Figure 1:**",
                "Barplot for pseudospectrum 27,",
                "displaying fragment *m/z* on the x-axis",
                "and intensity normalised to the maximum intensity",
                "on the y-axis.")
fig2 <- paste("**Figure 2:**",
                "Symmetric heat map of the distance matrix displaying",
                "pseudospectra similarities",
                "of the GC-EI-MS example data set",
                "along with dendrograms resulting from",
                "hierarchical clustering based on the distance matrix.",
                "The colour encoding is shown in the top-left insert.")
fig3 <- paste("**Figure 3:**",
                "Correlation network plot based on",
                "pseudospectra similarities",
                "of the GC-EI-MS example data set,",
                "generated with the default similarity threshold of `0.1`.",
                "Grey dots indicate non-identified features,",
                "orange dots identified ones.",
                "Labels display feature IDs, along with",
                "feature annotations, if existent.",
                "Edge widths are proportional to spectral similarity",
                "of the connected features.")
fig4 <- paste("**Figure 4:**",
                "Correlation network plot based on",
                "pseudospectra similarities",
                "of the GC-EI-MS example data set,",
                "generated with the custom similarity threshold of `0.4`.",
                "Grey dots indicate non-identified features,",
                "orange dots identified ones.",
                "Labels display feature IDs, along with",
                "feature annotations, if existent.",
                "Edge widths are proportional to spectral similarity",
                "of the connected features.")
fig5 <- paste("**Figure 5:**",
                "Circularised dendrogram as a result of",
                "agglomerative hierarchical clustering with average linkage",
                "as agglomeration criterion based on",
                "pseudospectra similarities",
                "of the GC-EI-MS example data set.",
                "Each leaf represents one feature and colours encode",
                "cluster affiliation of the features.",
                "Leaf labels display feature IDs, along with",
                "feature annotations, if existent.",
                "Distance from the central point is indicative",
                "of the height of the dendrogram.")
```

## Introduction

Although originally developed for high resolution LC-MS/MS data, 
CluMSID can also be used to find similarities in GC-EI-MS data, 
i.e. data from hard ionisation mass spectrometry.

As the peak picking and spectral merging differs considerably 
from data dependent ESI-MS/MS, we cannot use the standard `CluMSID` functions 
`extractMS2spectra()` and `mergeMS2spectra()`. 
In fact, the analysis of mass spectra from hard ionisation mass spectrometry 
resembles the one of MS^1^ pseudospectra in ESI-MS. 
Thus, we can use the CluMSID function `extractPseudospectra()` 
in conjunction with pseudspectra generated by the `CAMERA` package.

Since `xcms` and `CAMERA` sometimes have difficulties in handling 
GC-EI-MS data, we use the `metaMS` package that enables workflows 
specialised to the analysis of such data. We also require the `metaMSdata` 
package from which we import the `FEMSsettings` object that contains `xcms` 
and `CAMERA` settings for GC-EI-MS data.

```{r packages, message=FALSE, warning=FALSE}
library(CluMSID)
library(CluMSIDdata)
library(metaMS)
library(metaMSdata)
data(FEMsettings)
```

## Data import and preprocessing

As example data, we use GC-EI-MS metabolomics data from pooled cell extracts 
of *Pseudomonas aeruginosa* measured on a Thermo Scientific ITQ 
linear ion trap that has been converted to netCDF using Thermo Xcalibur. 
A netCDF file is available in the `CluMSIDdata` package:

```{r loading}
pool <- system.file("extdata", 
                    "1800802_TD_pool_total_1.cdf", 
                    package = "CluMSIDdata")
```

To generate a list of (pseudo)spectra, we first need an xsAnnotate object 
as generated by `CAMERA`. In the case of GC-MS data, it is more convenient 
to use to use the `metaMS` function `runCAMERA()` than actual `CAMERA` 
functions. `metaMS::runCAMERA` requires an `xcmsSet` object which we generate 
by using `xcms::xcmsSet` on our netCDF file (we can do that in one go). 
We used standard GC-MS settings for `runCAMERA()` as they are proposed in the 
[`metaMS` vignette](https://rdrr.io/bioc/metaMS/f/inst/doc/runGC.pdf).

```{r runCAMERA}
xA <- metaMS::runCAMERA(xcms::xcmsSet(pool), 
                chrom = "GC", 
                settings = metaMS::metaSetting(TSQXLS.GC, "CAMERA"))
```

## Extraction and annotation of spectra

From the `xsAnnotate` object, we can now extract the (pseudo)spectra 
using the `CluMSID` function `extractPseudospectra()` function 
as we would do for MS^1^ pseudospectra from LC-ESI-MS data.

```{r extract}
pslist <- extractPseudospectra(xA, min_peaks = 0)
```

Adding annotations is not as easy as with LC-(DDA-)MS/MS data, 
because only the retention time and the spectrum itself describe the feature 
and no precursor *m/z* is available. Thus, feature annotations/identifications 
made in a different programme, in this case MetaboliteDetector, 
have to be compared to the spectra in the `pslist` object.

Like with LC-(DDA-)MS/MS data, we can use `writeFeaturelist()` and 
`addAnnotations()` to add external annotations. 
The table output from `writeFeaturelist()` will give `NA` 
for all precursor *m/z*.

```{r Featurelist, eval=FALSE}
writeFeaturelist(pslist, "GC_pre.csv")
```

To facilitate manual annotation, it helps to plot the spectra 
along with the relevant information for every feature/pseudospectrum. 
That can be done by CluMSID's `specplot` function:

```{r temporary, fig.width=5, fig.asp=0.85, fig.cap=fig1}
specplot(pslist[[27]])
```

In this example, we load the list of feature annotations from `CluMSIDdata`:

```{r annotate}
apslist <- addAnnotations(featlist = pslist, 
                            annolist = system.file("extdata", 
                            "GC_post.csv", 
                            package = "CluMSIDdata"))
```

## Generation of distance matrix

This list of spectra in turn serves as an input for `distanceMatrix()`. 
As we are dealing with low resolution data, we have to adjust the 
*m/z* tolerance. The default value, 10ppm, is suitable for time-of-flight 
mass spectrometers while linear ion traps or single quadrupoles which are 
commonly used in GC-EI-MS only have unit mass resolution, equivalent to 
a relative mass error of 0.02 to 0.001 depending on the *m/z* of the analyte. 
We chose 0.02 to be tolerant enough for low molecular weight analytes:

```{r distmat}
pseudodistmat <- distanceMatrix(apslist, mz_tolerance = 0.02)
```

Starting from this distance matrix, we can use all the data exploration 
functions that `CluMSID` offers. In this example workflow, 
we look at a cluster dendrogram:

```{r heatmap, fig.width=7, fig.asp=1, fig.cap=fig2}
HCplot(pseudodistmat, type = "heatmap",
                cexRow = 0.4, cexCol = 0.4,
                margins = c(7,7))
```

It is directly visible that the resulting clusters are not as dense as with the
LC-MS/MS example data. In turn, there are more between-cluster similarities.
This also shows in the correlation network, resulting in a chaotic plot
when used with the default minimal similarity of `0.1`:

```{r network, fig.width=4.5, fig.asp=1, fig.cap=fig3}
networkplot(pseudodistmat, highlight_annotated = TRUE, 
                show_labels = TRUE, exclude_singletons = TRUE)
```

By choosing a higher similarity threshold of e.g. `0.4`, it is far easier
to identify clusters:

```{r network2, fig.width=4.5, fig.asp=1, fig.cap=fig4}
networkplot(pseudodistmat, highlight_annotated = TRUE, 
                show_labels = TRUE, exclude_singletons = TRUE,
                min_similarity = 0.4)
```

Presumably, the high between-cluster similarities are due to the low resolution
data and the resulting fact, that fragment with different chemical composition
but same unit resolution mass cannot be distinguished.

We can also use hierarchical clustering to identify clusters of similar
(pseudo-)spectra. Here, too, we have to adjust `h` to account for
higher between-cluster similarities:

```{r dendrogram, fig.width=7, fig.asp=1, fig.cap=fig5}
HCplot(pseudodistmat, h = 0.7, cex = 0.5)
```

We see that e.g. octadecanoic acid, hexadecanoic acid and dodecanoic acid
form a nice cluster as well as the phosphorate containing metabolites
phosphoenolpyruvic acid, glyceric acid-3-phosphate, glycerol-3-phosphate 
and phosphoric acid itself. It is also apparent that some features have
a similarity of 1 and could therefore represent the same compound, like
e.g. the features 98, 67 and 72. Those three features cluster together
with AMP and UMP, suggesting that they could be related to nucleotides.

To illustrate the use of CluMSID's accessory function with this
type of data, we take another look at nucleotides:
A signature fragment for nucleotides in GC-EI-MS is *m/z* 315 that derives
from pentose-5-phosphates. We see this fragment in Figure 1, the spectrum
of UMP (derivatised with 5 TMS groups). We can use findFragment to see
if there are more spectra outside the cluster that freature this fragment.
As we deal with unit masses, we would like to find *m/z* of 315 +/- 0.5
which we can do by setting `tolerance = 0.5/315`:

```{r fragment}
fragmentlist <- findFragment(apslist, mz = 315, tolerance = 0.5/315)

vapply(X = fragmentlist, FUN = accessID, FUN.VALUE = integer(1))
```

We find four more spectra that contain a 315 fragment that could be 
investigated closer.

In conclusion, every annotation method is extremely limited if only low 
resolution data is available and so is CluMSID. Still, we see that the tool
works independently of chromatography and mass spectrometry method
and even has the potential to give some good hints for feature annotation
in GC-EI-MS metabolomics.

# Session Info
```{r session}
sessionInfo()
```