synapter cross matching example

crossmatching.Rmd

# Cross Matching Example

```{r setup, include=FALSE}
library(knitr)
opts_chunk$set(fig.height=10, fig.width=14, fig.path="figure/crossmatching/", cache=TRUE, autodep=TRUE, tidy=FALSE)
```
## Default Synapter Preprocessing

(same as in `synergise`)

```{r loadsynapter}
suppressPackageStartupMessages(library("synapter"))
```

### File Import
```{r fileimport}
dataPath <- file.path("..", "..", "..", "playground", "data",
                      "biological_replicas_for_synapter")
iaPath <- file.path(dataPath, "HDMSE", "ia")
apexPath <- file.path(dataPath, "HDMSE", "apex")

input <- list(
  identpeptide=file.path(iaPath,
                         "fermentor_02_sample_15_HDMSE_01_20140303114821",
                         "fermentor_02_sample_15_HDMSE_01_IA_final_peptide.csv"),
  quantpeptide=file.path(iaPath,
                         "fermentor_03_sample_15_HDMSE_01_20140303120343",
                         "fermentor_03_sample_15_HDMSE_01_IA_final_peptide.csv"),
  quantpep3d=file.path(apexPath,
                       "fermentor_03_sample_15_HDMSE_01_20140302183822",
                       "fermentor_03_sample_15_HDMSE_01_Pep3DAMRT.csv"),
  fasta=file.path(dataPath,
                  "S.cerevisiae_uniprot_reference_canonical_16_09_13.fasta"))

## only for reproducibility
set.seed(1)
system.time(s <- Synapter(input))
```

### Filtering
```{r filtering}
filterUniqueDbPeptides(s, verbose=TRUE)
fdr <- fpr <- 0.05
ppm <- 20

filterIdentPepScore(s, method="BH", fdr=fdr)
filterQuantPepScore(s, method="BH", fdr=fdr)
filterIdentProtFpr(s, fpr=fpr)
filterQuantProtFpr(s, fpr=fpr)
filterIdentPpmError(s, ppm=ppm)
filterQuantPpmError(s, ppm=ppm)
```

### RT Model and Grid Search
```{r rtmodel}
mergePeptides(s)

modelRt(s, span=0.05)

system.time(searchGrid(s, verbose=FALSE))

setBestGridParams(s)

findEMRTs(s)
s
```

## Cross Matching

new workflow part

### File Import

#### Spectra
```{r spectraimport}
system.time(loadSpectrumXmlFiles(s, list(
  identspectrum=file.path(apexPath,
                          "fermentor_02_sample_15_HDMSE_01_20140302181938",
                          "fermentor_02_sample_15_HDMSE_01_Pep3D_Spectrum.xml"),
  quantspectrum=file.path(apexPath,
                          "fermentor_03_sample_15_HDMSE_01_20140302183822",
                          "fermentor_03_sample_15_HDMSE_01_Pep3D_Spectrum.xml")),
                        verbose=FALSE))
```

#### Fragments
```{r fragmentsimport}
system.time(loadFragmentCsvFiles(s, list(
  identfragments=file.path(iaPath,
                          "fermentor_02_sample_15_HDMSE_01_20140303114821",
                          "fermentor_02_sample_15_HDMSE_01_IA_final_fragment.csv"),
  quantfragments=file.path(iaPath,
                          "fermentor_03_sample_15_HDMSE_01_20140303120343",
                          "fermentor_03_sample_15_HDMSE_01_IA_final_fragment.csv")),
                          verbose=FALSE))
```

```{r importresults}
s
```

### Run Cross Matching
```{r cx}
## set mz tolerance (to consider two peaks as identical) to 25 ppm
setCrossMatchingPpmTolerance(s, 25)

crossMatching(s, verbose=FALSE)
```

### Plot Cross Matching Results
```{r cxsummary}
set.seed(1) # only for demonstration
crossMatchingSummaryMatrix <- plotCrossMatchingSummary(s)
set.seed(1) # only for demonstration
plotCrossMatchingSummary(s, matchColumn="spectrum.identXfragments.quant")
set.seed(1) # only for demonstration
plotCrossMatchingSummary(s, matchColumn="fragments.identXfragments.quant")
```
```{r cxsummarytable, results="asis"}
kable(crossMatchingSummaryMatrix, row.names=FALSE)
```

### Show Details Of Specific Data Points
find the `precursor.leID.ident` of the two false matches with more than 25
common peaks
```{r cxdetails}
cx <- s$CrossMatching
precursor.leID.ident <- cx$precursor.leID.ident[
 grepl("false", cx$matchType) & cx$fragments.identXfragments.quant > 25]

plotCrossMatching(s, key=precursor.leID.ident, column="precursor.leID.ident",
                  verbose=FALSE)
```
```{r cxtable, results="asis"}
options(digits=10)
kable(cx[cx$precursor.leID.ident %in% precursor.leID.ident,
      c("protein.Entry", "peptide.seq", "precursor.mhp", "precursor.retT",
        "precursor.leID.ident", "matched.quant.spectrumIDs",
        "fragments.identXfragments.quant")])
```
```{r cxdetails2}
matched.quant.spectrumIDs <-
  cx$matched.quant.spectrumIDs[cx$precursor.leID.ident %in% precursor.leID.ident]
pep3d <- s$QuantPep3DData
```
```{r pep3dtable, results="asis"}
kable(pep3d[pep3d$spectrumID %in% matched.quant.spectrumIDs,
      c("spectrumID", "mwHPlus", "ion_rt")])
```
```{r cxdetails3}
par(cex=5, mar=rep(0, 4))
plotFeatures(s, what="some", xlim=c(135, 137.7), ylim=c(1578.77, 1578.86))
plotFeatures(s, what="some", xlim=c(138.9, 141), ylim=c(1578.77, 1578.84))
```

Compare `"spectrum.identXfragments.quant` vs `"fragments.identXfragments.quant"`
and `"spectrum.quantXfragments.ident"` vs `"fragments.identXfragments.quant"`.

```{r cmpmatchcol}
plot(cx$fragments.identXfragments.quant, cx$spectrum.identXfragments.quant, col="#00000040", pch=19)
plot(cx$fragments.identXfragments.quant, cx$spectrum.quantXfragments.ident, col="#00000040", pch=19)
```

### Filtering Using The Cross Matching Results
`performance` before filtering
```{r precxfiltering}
performance(s)
```
cross match filtering
```{r cxfiltering}
setCrossMatchingMinimalNumberOfCommonPeaks(s, 3)
filterMatchedEMRTsByCommonPeaks(s)
```
`performance` after filtering
```{r postcxfiltering}
performance(s)
```