---
title: "Paired Mass Distance(PMD) analysis for GC/LC-MS based non-targeted analysis"
author: "Miao Yu"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{PMD analysis}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

## Introduction of Paired Mass Distance analysis

`pmd` package use Paired Mass Distance (PMD) relationship to analysis the GC/LC-MS based non-targeted data. PMD means the distance between two masses or mass to charge ratios. In mass spectrometry, PMD would keep the same between two masses or two mass to charge ratios(m/z). There are twe kinds of PMD involved in this package: PMD within same retention time group and PMD from different retention time groups.

### PMD within same retention time group

In GC/LC-MS based non-targeted analysis, peaks could be seperated by chromatograph. We could build retention time(RT) bins to assign peaks into different RT groups by retention time hierarchical clustering analysis. For each RT group, the peaks should come from same compounds or co-elutes. If certain PMD appeared in multiple RT groups, it would be related to the relationship about adducts, neutral loss, isotopologues or commen fragments ions.

### PMD from different retention time groups

The peaks from different retention time groups would like to be different compounds seperated by chromatograph. The PMD would reflect the relationship about homologous series or chemical reactions.

GlobalStd algorithm use the PMD within same RT group to find independent peaks among certain dataset. Structure/reaction directed analysis use PMD from different RT groups to screen important compounds or reactions.

## Data format

The input data should be a `list` object with at least two elements from a peaks list: 

- mass to charge ratio with name of `mz`, high resolution mass spectrometry is required
- retention time with name of `rt`

However, I suggested to add intensity and group information to the list for validation of PMD analysis.

In this package, a dataset from *in vivo* solid phase micro-extraction(SPME) was attached. This dataset contain 9 samples from 3 fish with triplicates samples for each fish. Here is the data strcture:

```{r demodata}
library(pmd)
data("spmeinvivo")
str(spmeinvivo)
```

You could build this `list` object from the `xcms` objects via `enviGCMS` package. When you have a `xcmsSet` object or `XCMSnExp` object named `xset`, you could use `enviGCMS::getmzrt(xset)` or `enviGCMS::getmzrt2(xset)` to get such list. Of course you could build such list by yourself.

## GlobalStd algorithm

GlobalStd algorithm try to find independent peaks among certain peaks list. The first step is retention time hierarchical clustering analysis. The second step is to find the relationship among adducts, neutral loss, isotopologues and commen fragments ions. The third step is to screen the independent peaks.

### Retention time hierarchical clustering

```{r rtg}
pmd <- getpaired(spmeinvivo, rtcutoff = 10, ng = 10)
plotrtg(pmd)
```

This plot would show the distribution of RT groups. The `rtcutoff` in function `getpaired` could be used to set the cutoff of the distances in retention time hierarchical clustering analysis.

### Relationship among adducts, neutral loss, isotopologues and commen fragments ions

The `ng` in function `getpaired` could be used to set cutoff of global PMD's retention time group numbers. If `ng` is 10, at least 10 of the retention time groups should contain the shown PMD relationship. You could use `plotpaired` to show the distribution.

```{r pmd}
plotpaired(pmd)
```

You could also show the distribution of PMD relationship by index:

```{r pmdindex}
# show the unique PMD found by getpaired function
for(i in 1:length(unique(pmd$paired$diff2))){
        diff <- unique(pmd$paired$diff2)[i]
        index <- pmd$paired$diff2 == diff
        plotpaired(pmd,index)
}
```

### Screen the independent peaks

You could use `getstd` function to get the independent peaks.

```{r std}
std <- getstd(pmd)
```

Here you could plot the peaks by `plotstd` function to show the distribution of independent peaks:

```{r stdplot}
plotstd(std)
```

You could also plot the peaks distribution by assign a retention time group via `plotstdrt`:

```{r stdrtplot}
par(mfrow = c(2,3))
plotstdrt(std,rtcluster = 23,main = 'Retention time group 23')
plotstdrt(std,rtcluster = 9,main = 'Retention time group 9')
plotstdrt(std,rtcluster = 18,main = 'Retention time group 18')
plotstdrt(std,rtcluster = 67,main = 'Retention time group 67')
plotstdrt(std,rtcluster = 49,main = 'Retention time group 49')
plotstdrt(std,rtcluster = 6,main = 'Retention time group 6')
```

### Use independent peaks for MS/MS validation

Independent peaks are supporsing generated from different compounds. We could use those peaks for MS/MS analysis instead of DIA or DDA. Here we need multiple injections for one sample since it might be impossible to get all ions' fragmental ions in one injection with good sensitivity. You could use `gettarget` to generate the index for the injections and output the peaks for each run.

```{r gettarget}
# you need retention time for independent peaks
index <- gettarget(std$rt[std$stdmassindex])
# output the ions for each injection
table(index)
# show the ions for the first injection
std$mz[index==1]
std$rt[index==1]
```


### Validation by principal components analysis(PCA)

You need to check the GlobalStd algorithm's results by principal components analysis(PCA).

```{r pca}
library(enviGCMS)
par(mfrow = c(1,2),mar = c(4,4,2,1)+0.1)
plotpca(std$data,lv = as.numeric(as.factor(std$group)),main = substitute(paste(italic('in vivo'), " SPME samples(all peaks)")))
plotpca(std$data[std$stdmassindex,],lv = as.numeric(as.factor(std$group)),main = substitute(paste(italic('in vivo'), " SPME samples(selected peaks)")))
```

### Comparision with other packages

GlobalStd algorithm in `pmd` package could be treated as a method to extract pseudospectra. You could use `getcluster` to get peaks groups information for all GlobalStd peaks. This function would consider the merge of GlobalStd peaks when certain peak is involved in multiple clusters. Then you could choose export peaks with the highest intensities in each GlobalStd merged peaks groups.

```{r comp}
stdcluster <- getcluster(std)
# extract pseudospectra for std peak 71
idx <- unique(stdcluster$cluster$largei[stdcluster$cluster$i==71])
plot(stdcluster$cluster$mz[stdcluster$cluster$largei==idx],stdcluster$cluster$ins[stdcluster$cluster$largei==idx],type = 'h',xlab = 'm/z',ylab = 'intensity',main = 'pseudospectra for GlobalStd peak 71')
# export peaks with the highest intensities in each GlobalStd peaks groups.
data <- stdcluster$data[stdcluster$stdmassindex2,]
```

You could also use `getcorcluster` to find peaks groups by correlation analysis only.

```{r corcomp}
corcluster <- getcorcluster(spmeinvivo)
par(mfrow = c(1,3),mar = c(4,4,2,1)+0.1)
plotpca(std$data,lv = as.numeric(as.factor(std$group)),main = substitute(paste(italic('in vivo'), " SPME samples(all peaks)")))
plotpca(std$data[std$stdmassindex,],lv = as.numeric(as.factor(std$group)),main = substitute(paste(italic('in vivo'), " SPME samples(selected peaks)")))
plotpca(std$data[corcluster$stdmassindex,],lv = as.numeric(as.factor(std$group)),main = substitute(paste(italic('in vivo'), " SPME samples(selected peaks by correlationship)")))
```


### GlobalStd algorithm with intensity data

GlobalStd algorithm is designed to analysis data without intensity data. However, if you have intensity data, the independant peaks could be selected with more confindence. You could set up cutoff of Pearson Correlation Coefficient between peaks to refine the peaks selected by GlobalStd within same retention time groups.

```{r globalcor}
std2 <- getstd(pmd,corcutoff = 0.9)
par(mfrow = c(1,3),mar = c(4,4,2,1)+0.1)
plotpca(std2$data,lv = as.numeric(as.factor(std2$group)),main = substitute(paste(italic('in vivo'), " SPME samples(all peaks)")))
plotpca(std$data[std$stdmassindex,],lv = as.numeric(as.factor(std$group)),main = substitute(paste(italic('in vivo'), " SPME samples(selected peaks)")))
plotpca(std2$data[std2$stdmassindex,],lv = as.numeric(as.factor(std2$group)),main = substitute(paste(italic('in vivo'), " SPME samples(selected peaks)")))
```


## Structure/Reaction directed analysis

`getsda` function could be used to perform Structure/reaction directed analysis. The cutoff of frequency is automately found by PMD network analysis with most cluster numbers.

```{r sda}
sda <- getsda(std)
```

You could use `plotstdsda` to show the distribution of the selected paired peaks.

```{r stdsda}
plotstdsda(sda)
```

You could also use index to show the distribution of certain PMDs.

```{r stdsdaidx}
par(mfrow = c(1,3),mar = c(4,4,2,1)+0.1)
plotstdsda(sda,sda$sda$diff2 == 2.02)
plotstdsda(sda,sda$sda$diff2 == 28.03)
plotstdsda(sda,sda$sda$diff2 == 58.04)
```

Structure/reaction directed analysis could be directily performed on all the peaks, which is slow to process:

```{r all}
sdaall <- getsda(spmeinvivo)
par(mfrow = c(1,3),mar = c(4,4,2,1)+0.1)
plotstdsda(sdaall,sdaall$sda$diff2 == 2.02)
plotstdsda(sdaall,sdaall$sda$diff2 == 28.03)
plotstdsda(sdaall,sdaall$sda$diff2 == 58.04)
```

### Structure/reaction directed analysis for peaks/compounds only

When you only have data of peaks without retention time or compounds list, structure/reaction directed analysis could also be done by `getrda` function.

```{r rda}
sda <- getrda(spmeinvivo$mz[std$stdmassindex])
```

### Structure/Reaction Network

One peak or compounds could be involved in multiple reactions. You could construct a network by such relationship.

If you have a specific compound and want to check the metabolites of certain PMD, you could use `getchain` to extract the network of that compounds

```{r tarnet}
library(igraph)
# check metabolites of C18H39NO
chain <- getchain(spmeinvivo,diff = c(2.02,14.02,15.99,58.04,13.98),mass = 286.3101,digits = 2,corcutoff = 0)
# show as network
net <- graph_from_data_frame(chain$sdac,directed = F)
pal <- (grDevices::colorRampPalette(rev(RColorBrewer::brewer.pal(11,"RdYlBu")
                )))(5)
plot(net,vertex.label=round(as.numeric(V(net)$name),2),vertex.size =5,edge.width = 5,edge.color = pal[as.numeric(as.factor(E(net)$diff2))],vertex.label.dist=1,vertex.color=ifelse(round(as.numeric(V(net)$name),4) %in% 286.3101,'red','black'),main = 'PMD network')
legend("topright",bty = "n",
       legend=unique(E(net)$diff2),
       fill=unique(pal[as.numeric(as.factor(E(net)$diff2))]), border=NA,horiz = F)
```

If you want to see all the independant peaks' high frequency PMDs as a network, the following code will help

```{r net}
sda <- getsda(std)
df <- sda$sda
net <- graph_from_data_frame(df,directed = F)
pal <- (grDevices::colorRampPalette(rev(RColorBrewer::brewer.pal(11,"RdYlBu")
                )))(length(unique(E(net)$diff2)))
plot(net,vertex.label=round(as.numeric(V(net)$name)),vertex.size = 7,edge.width = 5,edge.color = pal[as.numeric(as.factor(E(net)$diff2))],main = 'PMD network')
legend("topright",bty = "n",
       legend=unique(E(net)$diff2),
       fill=unique(pal[as.numeric(as.factor(E(net)$diff2))]), border=NA,horiz = F)
# Check the degree of the nodes
# Show the degree distribution of the vertices
deg <- degree(net, mode="all")
degree_distribution(net)
plot(net, vertex.size=deg/2,vertex.label=NA,vertex.size = 7, edge.width = 5)
# network community structure detection
ceb <- cluster_edge_betweenness(net,weights = abs(E(net)$cor), directed = F) 
plot(ceb, net,vertex.label=NA,) 
```

### Source appointment

Peaks from samples could be from endogenous compounds or exogenous compounds. However, it's hard to tell for untargeted analysis. In terms of PMD, if one peak belongs to a high frequency PMD network, it means a relatively high activity. If such sample belongs to a biological specimen, it might be endogenous compound. If a peak show no PMD network with other peaks, the biological system might not have enzyme to make reaction happen. Exogenous compounds will show a lower degree since they are xenbiotics. Since most of the peaks will show a low degree, the median of the degree could be used as cutoff. Then we could make source appointment if the assumption is hold.

```{r source}
median(deg)
endogenous <- names(deg)[deg>median(deg)]
exogenous <- names(deg)[deg<=median(deg)]
```

In this case, we will have `r length(endogenous)` endogenous compounds while `r length(exogenous)` exogenous compounds. When you find a peak show differences between groups, you could check the degree to infer its sources.

## Parameters selection

Retention time cluster cutoff should fit the peak picking algorithm. For HPLC, 10 is suggested and 5 could be used for UPLC.

Global PMD's retention time group numbers should be around 20 percent of the retention time cluster numbers. For example, if you find 100 retention time clusters, I suggested you use 20 as the empirical global PMD's retention time group numbers.

Another important hint is that pre-filter your peak list by black samples or other quality control samples. Otherwise the running time would be long and lots of pmd relationship would be just from noise.

## Wrap function

`globalstd` function is a wrap function to process GlobalStd algorithm and structure/reaction directed analysis in one line. All the plot function could be directly used on the `list` objects from `globalstd` function.

```{r wrap}
result <- globalstd(spmeinvivo,ng=10)
```

## Shiny application

An interactive document has been included in this package to perform PMD analysis. You need to prepare a csv file with m/z and retention time of peaks. Such csv file could be generated by run `enviGCMS::getmzrtcsv()` on the `list` object from `enviGCMS::getmzrt(xset)` or `enviGCMS::getmzrt2(xset)` function. You could also generate the csv file by `enviGCMS::getmzrt(xset,name = 'test')` or `enviGCMS::getmzrt2(xset, name = 'test')`. You will find the csv file in the working dictionary named `test.csv`.

Then you could run `runPMD()` to start the Graphical user interface(GUI) for GlobalStd algorithm and structure/reaction directed analysis. The PMD network for certain compound could be generated by run `runPMDnet()`.

## Database

To check the pmd reaction database:

```{r}
# all reaction
data("omics")
head(omics)
# kegg reaction
data("keggrall")
head(keggrall)
# literature reaction for mass spectrometry
data("sda")
head(sda)
```

To check the HMDB pmd database:

```{r}
data("hmdb")
head(hmdb)
```

## Conclusion

`pmd` package could be used to reduce the redundancy peaks for GC/LC-MS based research and perform structure/reaction directed analysis to screen known and unknown important compounds or reactions.

