2.1.1 Create a Synapter object

To demonstrate a typical step-by-step workflow we use example data that are available on http://proteome.sysbiol.cam.ac.uk/lgatto/synapter/data/. There is also an synobj2 object in synapterdata which contains the same data.

The Synapter constructor uses a named list of input files. Please note that we add identfragments (final_fragment.csv) and quantspectra (Spectrum.xml) because we want to apply the fragment matching later.

## Please find the raw data at:
## http://proteome.sysbiol.cam.ac.uk/lgatto/synapter/data/
library("synapter")

inlist <- list(
  identpeptide = "fermentor_03_sample_01_HDMSE_01_IA_final_peptide.csv.gz",
  identfragments = "fermentor_03_sample_01_HDMSE_01_IA_final_fragment.csv.gz",
  quantpeptide = "fermentor_02_sample_01_HDMSE_01_IA_final_peptide.csv.gz",
  quantpep3d = "fermentor_02_sample_01_HDMSE_01_Pep3DAMRT.csv.gz",
  quantspectra = "fermentor_02_sample_01_HDMSE_01_Pep3D_Spectrum.xml.gz",
  fasta = "S.cerevisiae_Uniprot_reference_canonical_18_03_14.fasta")

synobj2 <- Synapter(inlist, master=FALSE)

synobj2

## Object of class "Synapter" 
## Class version 2.0.0 
## Package version 1.99.0 
## Data files:
##  + Identification pep file: fermentor_03_sample_01_HDMSE_01_IA_final_peptide.csv.gz 
##  + Identification Fragment file: fermentor_02_sample_01_HDMSE_01_Pep3D_Spectrum.xml.gz 
##  + Quantitation pep file: fermentor_02_sample_01_HDMSE_01_IA_final_peptide.csv.gz 
##  + Quantitation Pep3DAMRT file: fermentor_02_sample_01_HDMSE_01_Pep3DAMRT.csv.gz 
##  + Quantitation Spectrum file: fermentor_02_sample_01_HDMSE_01_Pep3D_Spectrum.xml.gz 
##  + Fasta file: S.cerevisiae_Uniprot_reference_canonical_18_03_14.fasta 
## Log:
## [1] "Instance created on Tue Apr  4 09:01:34 2017"
## [2] "Read identification peptide data [11471,16]" 
## [ 10 lines ]
## [13] "Read identification fragment data [7078]"
## [14] "Read quantitation spectra [36539]"

2.1.2 Filtering

The first steps in each synapter analysis are filtering by peptide sequence, peptide length, ppm error and false positive rate.

Here we use the default values for each method. But the accompanying plotting methods should be used to find the best threshold:

filterUniqueDbPeptides(synobj2,
                       missedCleavages=0,
                       IisL=TRUE)
filterPeptideLength(synobj2, l=7)
plotFdr(synobj2)

filterQuantPepScore(synobj2, method="BH",
                    fdr=0.05)
filterIdentPepScore(synobj2, method="BH",
                    fdr=0.05)
par(mfcol=c(1, 2))
plotPpmError(synobj2, what="Ident")
plotPpmError(synobj2, what="Quant")

par(mfcol=c(1, 1))
filterQuantPpmError(synobj2, ppm=20)
filterIdentPpmError(synobj2, ppm=20)
plotPepScores(synobj2)

filterIdentProtFpr(synobj2, fpr=0.05)
filterQuantProtFpr(synobj2, fpr=0.05)

2.1.3 Modeling retention time

Next we merge the identified peptides from the identification run and quantification run and build a LOWESS based retention time model to remove systematic shifts in the retention times. Here we use the default values but as stated above the plotting methods should be used to find sensible thresholds.

mergePeptides(synobj2)

plotRt(synobj2, what="data")

setLowessSpan(synobj2, span=0.05)
modelRt(synobj2)

par(mfcol=c(1, 2))
plotRtDiffs(synobj2)
plotRt(synobj2, what="model", nsd=1)

par(mfcol=c(1, 1))

plotFeatures(synobj2, what="all", ionmobility=TRUE)

2.1.4 Grid search

To find EMRTS (exact m/z-retention time pairs) we try are running a grid search to find the best retention time tolerance and m/z tolerance that results in the most correct one-to-one matching in the merged (already identified) data. If the identification and quantitation run are HDMS\(^E\) data we could use the new 3D grid search that looks for the best matching in the retention time, m/z and ion mobility (drift time) domain to increase the accuracy. If one or both datasets are MS\(^E\) data it falls back to the traditional 2D grid search.

searchGrid(synobj2,
           imdiffs=seq(from=0.6, to=1.6, by=0.2),
           ppms=seq(from=2, to=20, by=2),
           nsds=seq(from=0.5, to=5, by=0.5))
setBestGridParams(synobj2)
findEMRTs(synobj2)
plotEMRTtable(synobj2)

2.1.5 Fragment matching

For the details of the fragment matching procedure we refer to the fragment matching vignette that is available online and with vignette("fragmentmatching", package = "synapter"). Briefly we compare the fragments of the identification run with the spectra from the quantification run and remove entries where there are very few/none common peaks/fragments between them.

First we starting by removing less intense fragments and peaks.

filterFragments(synobj2,
                what="fragments.ident",
                minIntensity=70)
filterFragments(synobj2,
                what="spectra.quant",
                minIntensity=70)

Next we look for common peaks via fragmentMatching:

fragmentMatching(synobj2)

## Warning in object$fragmentMatching(verbose = verbose): FragmentMatching ppm
## tolerance undefined. Setting to default value.

We get tables for unique and non-unique matches:

knitr::kable(fragmentMatchingPerformance(synobj2,
                                         what="unique"))

knitr::kable(fragmentMatchingPerformance(synobj2,
                                         what="non-unique"))

Subsequently we could filter by minimal accepted common peaks:

filterUniqueMatches(synobj2, minNumber=1)
filterNonUniqueMatches(synobj2, minDelta=2)
filterNonUniqueIdentMatches(synobj2)

Finally we rescue EMRTs that are filtered but were identified by PLGS:

rescueEMRTs(synobj2, method="rescue")

2.1.6 Modeling intensity

In a similar manner as correcting for the retention time drift we correct systematic errors of the intensity via a LOWESS model. The function modelIntensity has to applied after findEMRTs. The model is build on the merged peptides as it is done for the retention time model. But in contrast to the retention time model the prediction is necessary for the matched quantitation data.

plotIntensity(synobj2, what="data")

setLowessSpan(synobj2, 0.05)
modelIntensity(synobj2)
plotIntensity(synobj2, what="model", nsd=1)