examples.Rmd

```{r setup, cache=FALSE, include=FALSE}
library(knitr)
opts_chunk$set(fig.width=11, fig.height=7, base64_images=F, fig.align="center", message=F, warning=F, fig.path='figure_examples/')
```

# Example mz.unity usage
This example follows the analysis published in X.

### Load Data
```{r}
library(mz.unity)

spectrum = read.csv("glunad_spec.csv")

head(spectrum)
```
This spectrum was compiled from positive and negative mode scans from an LC/ESI/MS analysis of E. coli MG1655.  HILIC mode separation was employed with a Luna NH2 column coupled to a Q-Exactive Plus orbitrap mass spectrometer.

```{r}
plyr::count(spectrum$source)

ggplot(spectrum) + geom_segment(aes(x = mz.cont, xend = mz.cont, y = maxo, yend = 0)) + ggtitle("Spectrum from 21-22 minutes.")
```

### Format data for mz.unity

```{r}
mz = spectrum
mz$m = spectrum$mz
mz$z = sign(spectrum$mz.cont)
mz$d = 1
mz$source = spectrum$source
mz$id = seq(nrow(mz))
mz$row = seq(nrow(mz))

head(mz)
```

#### Find peaks with isotope support for Z = 2
```{r}
mz2 = mz; mz2$m = mz2$m * 2; mz2$z = mz2$z * 2
```

```{r}
ps = mz2[which(mz2$source %in% c("psp","psn")),]
pIso2 = mz.unity.search(A = ps, 
                        B = ps, 
                        M = M.iso, ppm = 1, 
                        BM.limits = cbind(M.min = c(1), M.max = c(1), B.n = c(1)))
pIso2 = reindex(pIso2, '^A$|^A.|^B$|^B.', ps$row)

head(pIso2)
```
Here we can see the z=2 isotope support detected in our spectrum.  In the first row, mz2[5,] plus granular formula M[1,] results in mz2[25,].  This is the M+1 carbon isotope of glutamate.

We do the same for background peaks:
```{r}

ps = mz2[which(mz2$source %in% c("bgp","bgn")),]
pIso2b = mz.unity.search(A = ps, 
                        B = ps, 
                        M = M.iso, ppm = 1, 
                        BM.limits = cbind(M.min = c(1), M.max = c(1), B.n = c(1)))
pIso2b = reindex(pIso2b, '^A$|^A.|^B$|^B.', ps$row)


mz2$monoiso = seq(nrow(mz2)) %in% c(pIso2$B.1, pIso2b$B.1) & !(seq(nrow(mz2)) %in% c(pIso2$A, pIso2b$A))
mz2$iso = seq(nrow(mz2)) %in% c(pIso2$B.1, pIso2b$B.1) | (seq(nrow(mz2)) %in% c(pIso2$A, pIso2b$A))
```
And record which ones have isotope support.  We assume all the rest are Z=1 peaks.

#### Find peaks with isotope support for Z = 1

```{r}

ps = mz[which(mz$source %in% c("psp","psn")),]
pIso = mz.unity.search(A = ps, 
                        B = ps, 
                        M = M.iso, ppm = 1, 
                        BM.limits = cbind(M.min = c(1), M.max = c(1), B.n = c(1)))
pIso = reindex(pIso, '^A$|^A.|^B$|^B.', ps$row)


ps = mz[which(mz$source %in% c("bgp","bgn")),]
pIsob = mz.unity.search(A = ps, 
                        B = ps, 
                        M = M.iso, ppm = 1, 
                        BM.limits = cbind(M.min = c(1), M.max = c(1), B.n = c(1)))
pIsob = reindex(pIsob, '^A$|^A.|^B$|^B.', ps$row)


mz$monoiso = seq(nrow(mz)) %in% c(pIso$B.1, pIsob$B.1) & !(seq(nrow(mz)) %in% c(pIso$A, pIsob$A))
mz$iso = seq(nrow(mz)) %in% c(pIso$B.1, pIsob$B.1) | (seq(nrow(mz)) %in% c(pIso$A, pIsob$A))

mzz = rbind(mz, subset(mz2, iso))
mzz$row = seq(mzz$row)
```


#### Annotate simple relationships

```{r}

ps = mzz[which(mzz$source %in% c("psp","psn") & ( mzz$monoiso | !mzz$iso )),]
pPol = mz.unity.search(A = ps, 
                       B = ps, 
                       M = M.z, ppm = 10, 
                       BM.limits = cbind(M.min = c(2), M.max = c(2), B.n = c(1))
                       )
pPol = reindex(pPol, '^M$|^M.', M.z$id)
pPol = subset(pPol, sign(ps[pPol$A, "z"]) != sign(ps[pPol$B.1,"z"]))
pPol = reindex(pPol, '^A$|^A.|^B$|^B.', ps$row)


ps = mzz[which(mzz$source %in% c("psp","psn") & ( mzz$monoiso | !mzz$iso )),]
pNl = mz.unity.search(A = ps, 
                      B = ps, 
                      M = M.n, ppm = 2, 
                      BM.limits = cbind(M.min = c(1), M.max = c(1), B.n = c(1))
                      )
pNl = reindex(pNl, '^M$|^M.', M.n$id)
pNl = reindex(pNl, '^A$|^A.|^B$|^B.', ps$row)


ps = mzz[which(mzz$source %in% c("psp","psn") & ( mzz$monoiso | !mzz$iso )),]
pZ = mz.unity.search(A = ps, 
                      B = ps, 
                      M = M.z, ppm = 2, 
                      BM.limits = cbind(M.min = c(1), M.max = c(2), B.n = c(1))
                      )
pZ = subset(pZ, sign(ps[pZ$A, "z"]) == sign(ps[pZ$B.1,"z"]))
pZ = reindex(pZ, '^A$|^A.|^B$|^B.', ps$row)

```



#### Annotate Analyte-Analyte Mers and Distal Fragments
```{r}
options(expressions=10000)

ps = mzz[which(mzz$source %in% c("psn") & ( mzz$monoiso | !mzz$iso )),] 
ms = subset(M.z, id == 11,drop=F)
pMern = mz.unity.search(A = ps, 
                        B = ps, 
                        M = ms, ppm = 2, 
                        BM.limits = cbind(M.min = c(1), M.max = c(1), B.n = c(2))
                        )
pMern = reindex(pMern, '^M$|^M.', ms$id)
pMern = reindex(pMern, '^A$|^A.|^B$|^B.', ps$row)


ps = mzz[which(mzz$source %in% c("psp") & ( mzz$monoiso | !mzz$iso )),] 
ms = subset(M.z, id == 11,drop=F)
pMerp = mz.unity.search(A = ps, 
                        B = ps, 
                        M = ms, ppm = 2, 
                        BM.limits = cbind(M.min = c(1), M.max = c(1), B.n = c(2))
                        )
pMerp = reindex(pMerp, '^M$|^M.', ms$id)
pMerp = reindex(pMerp, '^A$|^A.|^B$|^B.', ps$row)

pMer = rbind(pMerp, pMern)
```

#### Considering only peaks with isotope support - Annotate Analyte-Analyte Mers and Distal Fragments

```{r}
options(expressions=10000)

ps = mzz[which(mzz$source %in% c("psn") & ( mzz$monoiso )),] 
ms = subset(M.z, id == 11,drop=F)
pMern = mz.unity.search(A = ps, 
                        B = ps, 
                        M = ms, ppm = 2, 
                        BM.limits = cbind(M.min = c(1), M.max = c(1), B.n = c(2))
                        )
pMern = reindex(pMern, '^M$|^M.', ms$id)
pMern = reindex(pMern, '^A$|^A.|^B$|^B.', ps$row)


ps = mzz[which(mzz$source %in% c("psp") & ( mzz$monoiso )),] 
ms = subset(M.z, id == 11,drop=F)
pMerp = mz.unity.search(A = ps, 
                        B = ps, 
                        M = ms, ppm = 2, 
                        BM.limits = cbind(M.min = c(1), M.max = c(1), B.n = c(2))
                        )
pMerp = reindex(pMerp, '^M$|^M.', ms$id)
pMerp = reindex(pMerp, '^A$|^A.|^B$|^B.', ps$row)

pMerI = rbind(pMerp, pMern)
```


##### Annotate Cross Polarity Mers and Distal Fragments
```{r}
options(expressions=10000)

ps = mzz[which(mzz$source %in% c("psp", "psn") & ( mzz$monoiso | !mzz$iso )),] 
ms = subset(M.z, id == 11,drop=F)
pFrag = mz.unity.search(A = ps, 
                        B = ps, 
                        M = ms, ppm = 2, 
                        BM.limits = cbind(M.min = c(1), M.max = c(1), B.n = c(2))
                        )
pFrag = subset(pFrag, xor(ps$source[B.1] == "psp", ps$source[B.2] == "psp"))
pFrag = reindex(pFrag, '^M$|^M.', ms$id)
pFrag = reindex(pFrag, '^A$|^A.|^B$|^B.', ps$row)
```


####  Considering only peaks with isotope support - Annotate Cross Polarity Mers and Distal Fragments
```{r}
options(expressions=10000)

ps = mzz[which(mzz$source %in% c("psp", "psn") & ( mzz$monoiso )),] 
ms = subset(M.z, id == 11,drop=F)
pFragI = mz.unity.search(A = ps, 
                        B = ps, 
                        M = ms, ppm = 2, 
                        BM.limits = cbind(M.min = c(1), M.max = c(1), B.n = c(2))
                        )
pFragI = subset(pFragI, xor(ps$source[B.1] == "psp", ps$source[B.2] == "psp"))
pFragI = reindex(pFragI, '^M$|^M.', ms$id)
pFragI = reindex(pFragI, '^A$|^A.|^B$|^B.', ps$row)
```

##### Annotate Background Mers
```{r}
options(expressions=40000)

psa = mzz[which(mzz$source %in% c("psn") & ( mzz$monoiso | !mzz$iso )),] 
psb = rbind(
  mzz[which(mzz$source %in% c("psn") & ( mzz$monoiso | !mzz$iso )),],
  mzz[which(mzz$source %in% c("bgn") & ( mzz$monoiso )),]
  )
ms = subset(M.z, id == 11,drop=F)
pMernb = mz.unity.search(A = psa, 
                        B = psb, 
                        M = ms, ppm = 2, 
                        BM.limits = cbind(M.min = c(1), M.max = c(1), B.n = c(2))
                        )
pMernb = subset(pMernb, xor(psb$source[B.1] == "psn", psb$source[B.2] == "psn"))
pMernb = reindex(pMernb, '^M$|^M.', ms$id)
pMernb = reindex(pMernb, '^A$|^A.', psa$row)
pMernb = reindex(pMernb, '^B$|^B.', psb$row)


psa = mzz[which(mzz$source %in% c("psp") & ( mzz$monoiso | !mzz$iso )),] 
psb = rbind(
  mzz[which(mzz$source %in% c("psp") & ( mzz$monoiso | !mzz$iso )),],
  mzz[which(mzz$source %in% c("bgp") & ( mzz$monoiso )),]
  )
ms = subset(M.z, id == 11,drop=F)
pMerpb = mz.unity.search(A = psa, 
                        B = psb, 
                        M = ms, ppm = 2, 
                        BM.limits = cbind(M.min = c(1), M.max = c(1), B.n = c(2))
                        )
pMerpb = subset(pMerpb, xor(psb$source[B.1] == "psp", psb$source[B.2] == "psp"))
pMerpb = reindex(pMerpb, '^M$|^M.', ms$id)
pMerpb = reindex(pMerpb, '^A$|^A.', psa$row)
pMerpb = reindex(pMerpb, '^B$|^B.', psb$row)

pBgMer = rbind(pMernb, pMerpb)
```
#### Considering only peaks with isotope support - Annotate Background Mers
```{r}
options(expressions=20000)

psa = mzz[which(mzz$source %in% c("psn") & ( mzz$monoiso )),] 
psb = rbind(
  mzz[which(mzz$source %in% c("psn") & ( mzz$monoiso )),],
  mzz[which(mzz$source %in% c("bgn") & ( mzz$monoiso )),]
  )
ms = subset(M.z, id == 11,drop=F)
pMernb = mz.unity.search(A = psa, 
                        B = psb, 
                        M = ms, ppm = 2, 
                        BM.limits = cbind(M.min = c(1), M.max = c(1), B.n = c(2))
                        )
pMernb = subset(pMernb, xor(psb$source[B.1] == "psn", psb$source[B.2] == "psn"))
pMernb = reindex(pMernb, '^M$|^M.', ms$id)
pMernb = reindex(pMernb, '^A$|^A.', psa$row)
pMernb = reindex(pMernb, '^B$|^B.', psb$row)


psa = mzz[which(mzz$source %in% c("psp") & ( mzz$monoiso )),] 
psb = rbind(
  mzz[which(mzz$source %in% c("psp") & ( mzz$monoiso )),],
  mzz[which(mzz$source %in% c("bgp") & ( mzz$monoiso )),]
  )
ms = subset(M.z, id == 11,drop=F)
pMerpb = mz.unity.search(A = psa, 
                        B = psb, 
                        M = ms, ppm = 2, 
                        BM.limits = cbind(M.min = c(1), M.max = c(1), B.n = c(2))
                        )
pMerpb = subset(pMerpb, xor(psb$source[B.1] == "psp", psb$source[B.2] == "psp"))
pMerpb = reindex(pMerpb, '^M$|^M.', ms$id)
pMerpb = reindex(pMerpb, '^A$|^A.', psa$row)
pMerpb = reindex(pMerpb, '^B$|^B.', psb$row)

pBgMerI = rbind(pMernb, pMerpb)
```

### Aggregate Relationships
```{r}
pZ$rel = "z"
pPol$rel = "pol"
pNl$rel = "nl"
pMer$rel = "mer"
pMerI$rel = "meri"
pFrag$rel = "frag"
pFragI$rel = "fragi"
pBgMer$rel = "bg"
pBgMerI$rel = "bgi"

rels = do.call(plyr::rbind.fill, list(
  pZ,
  pPol,
  pNl,
  pMer,
  pFrag,
  pBgMer,
  pBgMerI,
  pFragI,
  pMerI
  ))

```

```{r}
plyr::count(mz$source)
```
Peak relationships with Z = 1
```{r}
plyr::count(mzz$source)
```
Peak relationships with Z = 1 and 2

```{r}
plyr::count(rels$rel)
```
Types of relationships detected

## Plot Charge Carriers

```{r}

df = subset(rels, rel %in% c("z", "pol"))
df = reindex(df, '^A$|^A.|^B$|^B.', mzz$id)

vertices = c(unique(subset(mzz, monoiso & source %in% c("psn", "psp"))[,"id"]))
vertices = unique(c(vertices, unlist(expandGraph(df)[,1:2])))

g = graph.data.frame(expandGraph(df)[,1:2], vertices = data.frame(v=as.character(vertices)))
plot(g, asp = .4, vertex.size = 2, edge.arrow.mode=1, edge.arrow.size = 0.3, vertex.label.cex=.6, vertex.frame.color='transparent', vertex.label.color = "transparent")

length(decompose(g))
```


## Aggregate and Plot Fine Isotopic Spectrum from relationship graph
```{r}
isos = aggregate.children(pIso, 1)

df = mzz[isos,]
df$group = round(df$mz)
ggplot(df) + geom_segment(aes(x = mz, xend = mz, y = maxo, yend = 0)) + facet_grid(.~group, scales="free")
```

## Examine ambiguous neutral losses

```{r}
df = subset(rels, rel %in% "nl")
df$ambig = df$M.1 %in% c(16,18) 
nrow(df)

cnts = plyr::count(df$M.1)
cnts$name = rownames(M.n)[match(cnts$x,M.n$id)]
cnts

plyr::count(subset(df, rel %in% "nl", "ambig"))
```

## Plot all simple relationships
```{r}

df = subset(rels, rel %in% c("z", "pol", "nl"))
df = reindex(df, '^A$|^A.|^B$|^B.', mzz$id)

vertices = c(unique(subset(mzz, monoiso & source %in% c("psn", "psp"))[,"id"]))
vertices = unique(c(vertices, unlist(expandGraph(df)[,1:2])))

g = graph.data.frame(expandGraph(df)[,1:2], vertices = data.frame(v=as.character(vertices)))
plot(g, asp = .4, vertex.size = 2, edge.arrow.mode=1, edge.arrow.size = 0.3, vertex.label.cex=.6, vertex.frame.color='transparent', vertex.label.color = "transparent")

edgetypes = cbind(
  rel = c("nl", "pol", "z"),
  label = c("NL", "Polarity", "Charge"),
  color = c("cornflowerblue", "red", "darkgoldenrod1")
  )

spec = mzz[,c("mz.cont", "id", "maxo")]
colnames(spec) =c("mz", "id", "int")

dft = cbind(df, g = aggregate.self(df, F))

plyr::count(dft$g.A)

subset(dft, A %in% c(1, 5, 9))
subset(dft, B.1 %in% c(1, 5, 9))

dft2 = subset(dft, g.A %in% c(1, 2, 4))
unity.graph(dft2, spec, edgetypes)

```


## Plot of relationships including mers and fragments
Note that not all of these relationships imply that all peaks involved are from the same compound.
```{r}

df = subset(rels, rel %in% c("z", "pol", "nl", "frag", "mer"))
df = reindex(df, '^A$|^A.|^B$|^B.', mzz$id)

vertices = c(unique(subset(mzz, monoiso & source %in% c("psn", "psp"))[,"id"]), 22)
vertices = unique(c(vertices, unlist(expandGraph(df)[,1:2])))

g = graph.data.frame(expandGraph(df)[,1:2], vertices = data.frame(v=as.character(vertices)))
plot(g, asp = .4, vertex.size = 2, edge.arrow.mode=1, edge.arrow.size = 0.3, vertex.label.cex=.6, vertex.frame.color='transparent', vertex.label.color = "transparent")

length(decompose(g))
```


## Explore Glutamate and NAD Cluster
```{r}
tmprels1 = subset(rels, rel %in% c("nl", "pol", "z"))
tmprels2 = subset(rels, rel %in% c("xmer", "frag"))
tmprels = rbind(tmprels1, tmprels2)
tmprels = cbind(tmprels, g = aggregate.self(tmprels))
head(tmprels)
```
Columns starting with g are the compound groups.  This assumes that all detected mers are mers and not fragments and that no ambiguous relationships connect distinct analytes.


```{r}
gi = matrix(c(unlist(tmprels[,grep('^A|^B.', colnames(tmprels))]), unlist(tmprels[,grep('^g.', colnames(tmprels))])), ncol = 2, dimnames = list(NULL, c("i", "g")))
gi = unique(data.frame(gi))
gi = gi[!is.na(gi$g),]

g.glu = subset(gi, i == 1)$g
g.nad = subset(gi, i == 5)$g
g.glunad = subset(gi, i == 9)$g
gs = c(g.glu, g.nad, g.glunad)

head(gi)
```
There are many groups but we are only going to plot those containing glutamate, nad and their dimer.

```{r}
agroup = subset(gi, g %in% gs)
rels.1 = subset(tmprels, (g.A %in% gs & g.B.1 %in% gs & (g.B.2 %in% gs | is.na(g.B.2))))
rels.1$rel[rels.1$rel == "mer" & rels.1$g.A == rels.1$g.B.1] = "mer.self"

es = expandGraph(rels.1)
es$rel = rels.1[es$row, "rel"]
es$color = c("lightgreen", "lightgreen", "orange", "blueviolet", "blueviolet", "blue", "red", "blue", "blue", "blue")[as.numeric(factor(es$rel))]
vs = unique(unlist(es[,1:2]))

es.l = expandGraph(subset(rels.1, g.A == g.B.1))
g.l = graph.data.frame(es.l, vertices = vs)
l = layout.auto(g.l)

g = graph.data.frame(es, vertices = vs)
plot(g, asp = .4, vertex.size = 2, edge.arrow.mode=1, edge.arrow.size = 0.3, vertex.label.cex=.6, vertex.frame.color='transparent', vertex.label.color = "transparent", layout = l)


mat = match(as.numeric(V(g)$name), gi$i)
c2 = list(membership = gi$g[mat], names = as.character(gi$i)[mat])
class(c2) = "communities"
plot(c2, g, asp = .4, vertex.size = 2, edge.arrow.mode=1, edge.arrow.size = 0.3, vertex.label.cex=.6, vertex.frame.color='transparent', vertex.label.color = "transparent", layout = l)
```