Mass spectrometry and proteomics data analysis

1 Setup

The follow packages will be used throughout this documents. R version 3.3.1 or higher is required to install all the packages using BiocInstaller::biocLite.

library("mzR")
library("mzID")
library("MSnID")
library("MSnbase")
library("rpx")
library("MLInterfaces")
library("pRoloc")
library("pRolocdata")
library("MSGFplus")
library("rols")
library("hpar")

The most convenient way to install all the tutorials requirement (and more related content), is to install RforProteomics with all its dependencies.

library("BiocInstaller")
biocLite("RforProteomics", dependencies = TRUE)

Other packages of interest, such as rTANDEM or MSGFgui will be described later in the document but are not required to execute the code in this workflow.

2 Introduction

This workflow illustrates R / Bioconductor infrastructure for proteomics. Topics covered focus on support for open community-driven formats for raw data and identification results, packages for peptide-spectrum matching, data processing and analysis:

  • Exploring available infrastructure
  • Mass spectrometry data
  • Getting data from proteomics repositories
  • Handling raw MS data
  • Handling identification data
  • MS/MS database search
  • Analysing search results
  • High-level data interface
  • Quantitative proteomics
  • Importing third-party quantitation data
  • Data processing and analysis
  • Statistical analysis
  • Machine learning
  • Annotation
  • Other relevant packages/pipelines

Links to other packages and references are also documented. In particular, the vignettes included in the RforProteomics package also contains relevant material.

3 Exploring available infrastructure

In Bioconductor version 3.5, there are respectively 94 proteomics, 64 mass spectrometry software packages and 17 mass spectrometry experiment packages. These respective packages can be extracted with the proteomicsPackages(), massSpectrometryPackages() and massSpectrometryDataPackages() and explored interactively.

library("RforProteomics")
pp <- proteomicsPackages()
display(pp)

4 Mass spectrometry data

Most community-driven formats are supported in R, as detailed in the table below.

Type Format Package
raw mzML, mzXML, netCDF, mzData mzR (read)
identification mzIdentML mzR (read) and mzID (read)
quantitation mzQuantML
peak lists mgf MSnbase (read/write)
other mzTab MSnbase (read)

5 Getting data from proteomics repositories

MS-based proteomics data is disseminated through the ProteomeXchange infrastructure, which centrally coordinates submission, storage and dissemination through multiple data repositories, such as the PRIDE data base at the EBI for MS/MS experiments, PASSEL at the ISB for SRM data and the MassIVE resource. The rpx is an interface to ProteomeXchange and provides a basic access to PX data.

library("rpx")
pxannounced()
## 15 new ProteomeXchange annoucements
##     Data.Set    Publication.Data             Message
## 1  PXD007078 2017-07-21 09:00:55                 New
## 2  PXD005852 2017-07-21 08:40:16 Updated information
## 3  PXD005139 2017-07-21 08:29:09 Updated information
## 4  PXD005140 2017-07-21 08:28:42 Updated information
## 5  PXD006754 2017-07-21 07:43:03                 New
## 6  PXD006356 2017-07-21 07:39:22 Updated information
## 7  PXD006925 2017-07-21 07:38:50                 New
## 8  PXD003360 2017-07-21 07:35:45                 New
## 9  PXD005529 2017-07-20 13:49:05                 New
## 10 PXD006598 2017-07-20 11:36:25                 New
## 11 PXD005139 2017-07-20 11:14:41                 New
## 12 PXD005140 2017-07-20 11:12:39                 New
## 13 PXD006356 2017-07-20 11:07:15                 New
## 14 PXD002295 2017-07-20 10:57:32 Updated information
## 15 PXD002673 2017-07-20 10:53:45 Updated information

Using the unique PXD000001 identifier, we can retrieve the relevant metadata that will be stored in a PXDataset object. The names of the files available in this data can be retrieved with the pxfiles accessor function.

px <- PXDataset("PXD000001")
px
## Object of class "PXDataset"
##  Id: PXD000001 with 12 files
##  [1] 'F063721.dat' ... [12] 'generated'
##  Use 'pxfiles(.)' to see all files.
pxfiles(px)
##  [1] "F063721.dat"                                                         
##  [2] "F063721.dat-mztab.txt"                                               
##  [3] "PRIDE_Exp_Complete_Ac_22134.xml.gz"                                  
##  [4] "PRIDE_Exp_mzData_Ac_22134.xml.gz"                                    
##  [5] "PXD000001_mztab.txt"                                                 
##  [6] "README.txt"                                                          
##  [7] "TMT_Erwinia_1uLSike_Top10HCD_isol2_45stepped_60min_01-20141210.mzML" 
##  [8] "TMT_Erwinia_1uLSike_Top10HCD_isol2_45stepped_60min_01-20141210.mzXML"
##  [9] "TMT_Erwinia_1uLSike_Top10HCD_isol2_45stepped_60min_01.mzXML"         
## [10] "TMT_Erwinia_1uLSike_Top10HCD_isol2_45stepped_60min_01.raw"           
## [11] "erwinia_carotovora.fasta"                                            
## [12] "generated"

Other metadata for the px data set:

pxtax(px)
## [1] "Erwinia carotovora"
pxurl(px)
## [1] "ftp://ftp.pride.ebi.ac.uk/pride/data/archive/2012/03/PXD000001"
pxref(px)
## [1] "<!DOCTYPE HTML PUBLIC \"-//IETF//DTD HTML 2.0//EN\"> <html><head> <title>301 Moved Permanently</title> </head><body> <h1>Moved Permanently</h1> <p>The document has moved <a href=\"https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi/efetch.fcgi?db=pubmed&amp;id=23692960&amp;rettype=docsum&amp;retmode=text\">here</a>.</p> </body></html>"

Data files can then be downloaded with the pxget function. Below, we retrieve the raw data file. The file is downloaded in the working directory and the name of the file is return by the function and stored in the mzf variable for later use.

fn <- "TMT_Erwinia_1uLSike_Top10HCD_isol2_45stepped_60min_01-20141210.mzML"
mzf <- pxget(px, fn)
## Downloading 1 file
mzf
## [1] "/liftrroot/TMT_Erwinia_1uLSike_Top10HCD_isol2_45stepped_60min_01-20141210.mzML"

6 Handling raw MS data

The mzR package provides an interface to the proteowizard C/C++ code base to access various raw data files, such as mzML, mzXML, netCDF, and mzData. The data is accessed on-disk, i.e it is not loaded entirely in memory by default but only when explicitly requested. The three main functions are openMSfile to create a file handle to a raw data file, header to extract metadata about the spectra contained in the file and peaks to extract one or multiple spectra of interest. Other functions such as instrumentInfo, or runInfo can be used to gather general information about a run.

Below, we access the raw data file downloaded in the previous section and open a file handle that will allow us to extract data and metadata of interest.

library("mzR")
ms <- openMSfile(mzf)
ms
## Mass Spectrometry file handle.
## Filename:  TMT_Erwinia_1uLSike_Top10HCD_isol2_45stepped_60min_01-20141210.mzML 
## Number of scans:  7534

The header function returns the metadata of all available peaks:

hd <- header(ms)
dim(hd)
## [1] 7534   21
names(hd)
##  [1] "seqNum"                   "acquisitionNum"          
##  [3] "msLevel"                  "polarity"                
##  [5] "peaksCount"               "totIonCurrent"           
##  [7] "retentionTime"            "basePeakMZ"              
##  [9] "basePeakIntensity"        "collisionEnergy"         
## [11] "ionisationEnergy"         "lowMZ"                   
## [13] "highMZ"                   "precursorScanNum"        
## [15] "precursorMZ"              "precursorCharge"         
## [17] "precursorIntensity"       "mergedScan"              
## [19] "mergedResultScanNum"      "mergedResultStartScanNum"
## [21] "mergedResultEndScanNum"

We can extract metadata and scan data for scan 1000 as follows:

hd[1000, ]
##      seqNum acquisitionNum msLevel polarity peaksCount totIonCurrent
## 1000   1000           1000       2        0        274       1048554
##      retentionTime basePeakMZ basePeakIntensity collisionEnergy
## 1000      1106.916    136.061            164464              45
##      ionisationEnergy    lowMZ   highMZ precursorScanNum precursorMZ
## 1000                0 104.5467 1370.758              992    683.0817
##      precursorCharge precursorIntensity mergedScan mergedResultScanNum
## 1000               2           689443.7          0                   0
##      mergedResultStartScanNum mergedResultEndScanNum
## 1000                        0                      0
head(peaks(ms, 1000))
##          [,1]     [,2]
## [1,] 104.5467 308.9326
## [2,] 104.5684 308.6961
## [3,] 108.8340 346.7183
## [4,] 109.3928 365.1236
## [5,] 110.0345 616.7905
## [6,] 110.0703 429.1975
plot(peaks(ms, 1000), type = "h")

Below we reproduce the example from the MSmap function from the MSnbase package to plot a specific slice of the raw data using the mzR functions we have just described.

## a set of spectra of interest: MS1 spectra eluted
## between 30 and 35 minutes retention time
ms1 <- which(hd$msLevel == 1)
rtsel <- hd$retentionTime[ms1] / 60 > 30 &
    hd$retentionTime[ms1] / 60 < 35

## the map
M <- MSmap(ms, ms1[rtsel], 521, 523, .005, hd)
## 1
plot(M, aspect = 1, allTicks = FALSE)

plot3D(M)

## With some MS2 spectra
i <- ms1[which(rtsel)][1]
j <- ms1[which(rtsel)][2]
M2 <- MSmap(ms, i:j, 100, 1000, 1, hd)
## 1
plot3D(M2)

7 Handling identification data

The RforProteomics package distributes a small identification result file (see ?TMT_Erwinia_1uLSike_Top10HCD_isol2_45stepped_60min_01.mzid) that we load and parse using infrastructure from the mzID package.

library("mzID")
f <- dir(system.file("extdata", package = "RforProteomics"),
         pattern = "mzid", full.names=TRUE)
basename(f)
## [1] "TMT_Erwinia.mzid.gz"
id <- mzID(f)
## reading TMT_Erwinia.mzid.gz... DONE!
id
## An mzID object
## 
## Software used:   MS-GF+ (version: Beta (v10072))
## 
## Rawfile:         /home/lgatto/dev/00_github/RforProteomics/sandbox/TMT_Erwinia_1uLSike_Top10HCD_isol2_45stepped_60min_01.mzXML
## 
## Database:        /home/lgatto/dev/00_github/RforProteomics/sandbox/erwinia_carotovora.fasta
## 
## Number of scans: 5287
## Number of PSM's: 5563

Various data can be extracted from the mzID object, using one the accessor functions such as database, scans, peptides, … The object can also be converted into a data.frame using the flatten function.

The mzR package also provides support fasta parsing mzIdentML files with the openIDfile function. As for raw data, the underlying C/C++ code comes from the proteowizard.

library("mzR")
f <- dir(system.file("extdata", package = "RforProteomics"),
         pattern = "mzid", full.names=TRUE)

id1 <- openIDfile(f)
fid1 <- mzR::psms(id1)

head(fid1)
##   spectrumID chargeState rank passThreshold experimentalMassToCharge
## 1  scan=5782           3    1          TRUE                1080.2325
## 2  scan=6037           3    1          TRUE                1002.2089
## 3  scan=5235           3    1          TRUE                1189.2836
## 4  scan=5397           3    1          TRUE                 960.5365
## 5  scan=6075           3    1          TRUE                1264.3409
## 6  scan=5761           2    1          TRUE                1268.6429
##   calculatedMassToCharge                            sequence modNum
## 1              1080.2321 PVQIQAGEDSNVIGALGGAVLGGFLGNTIGGGSGR      0
## 2              1002.2115        TQVLDGLINANDIEVPVALIDGEIDVLR      0
## 3              1189.2800   TKGLNVMQNLLTAHPDVQAVFAQNDEMALGALR      0
## 4               960.5365         SQILQQAGTSVLSQANQVPQTVLSLLR      0
## 5              1264.3419 PIIGDNPFVVVLPDVVLDESTADQTQENLALLISR      0
## 6              1268.6501             WTSQSSLDLGEPLSLITESVFAR      0
##   isDecoy post pre start end DatabaseAccess DBseqLength DatabaseSeq
## 1   FALSE    S   R    50  84        ECA1932         155            
## 2   FALSE    R   K   288 315        ECA1147         434            
## 3   FALSE    A   R   192 224        ECA0013         295            
## 4   FALSE    -   R   264 290        ECA1731         290            
## 5   FALSE    F   R   119 153        ECA1443         298            
## 6   FALSE    Y   K   264 286        ECA1444         468            
##                                         DatabaseDescription acquisitionNum
## 1                        ECA1932 outer membrane lipoprotein           5782
## 2                                    ECA1147 trigger factor           6037
## 3                ECA0013 ribose-binding periplasmic protein           5235
## 4                                         ECA1731 flagellin           5397
## 5      ECA1443 UTP--glucose-1-phosphate uridylyltransferase           6075
## 6 ECA1444 6-phosphogluconate dehydrogenase, decarboxylating           5761

9 Analysing search results

The MSnID package can be used for post-search filtering of MS/MS identifications. One starts with the construction of an MSnID object that is populated with identification results that can be imported from a data.frame or from mzIdenML files. Here, we will use the example identification data provided with the package.

mzids <- system.file("extdata", "c_elegans.mzid.gz", package="MSnID")
basename(mzids)
## [1] "c_elegans.mzid.gz"

We start by loading the package, initialising the MSnID object, and add the identification result from our mzid file (there could of course be more that one).

library("MSnID")
msnid <- MSnID(".")
## Note, the anticipated/suggested columns in the
## peptide-to-spectrum matching results are:
## -----------------------------------------------
## accession
## calculatedMassToCharge
## chargeState
## experimentalMassToCharge
## isDecoy
## peptide
## spectrumFile
## spectrumID
msnid <- read_mzIDs(msnid, mzids)
## Reading from mzIdentMLs ...
## reading c_elegans.mzid.gz... DONE!
show(msnid)
## MSnID object
## Working directory: "."
## #Spectrum Files:  1 
## #PSMs: 12263 at 36 % FDR
## #peptides: 9489 at 44 % FDR
## #accessions: 7414 at 76 % FDR

Printing the MSnID object returns some basic information such as

  • Working directory.
  • Number of spectrum files used to generate data.
  • Number of peptide-to-spectrum matches and corresponding FDR.
  • Number of unique peptide sequences and corresponding FDR.
  • Number of unique proteins or amino acid sequence accessions and corresponding FDR.

The package then enables to define, optimise and apply filtering based for example on missed cleavages, identification scores, precursor mass errors, etc. and assess PSM, peptide and protein FDR levels. To properly function, it expects to have access to the following data

## [1] "accession"                "calculatedMassToCharge"  
## [3] "chargeState"              "experimentalMassToCharge"
## [5] "isDecoy"                  "peptide"                 
## [7] "spectrumFile"             "spectrumID"

which are indeed present in our data:

names(msnid)
##  [1] "spectrumID"                "scan number(s)"           
##  [3] "acquisitionNum"            "passThreshold"            
##  [5] "rank"                      "calculatedMassToCharge"   
##  [7] "experimentalMassToCharge"  "chargeState"              
##  [9] "MS-GF:DeNovoScore"         "MS-GF:EValue"             
## [11] "MS-GF:PepQValue"           "MS-GF:QValue"             
## [13] "MS-GF:RawScore"            "MS-GF:SpecEValue"         
## [15] "AssumedDissociationMethod" "IsotopeError"             
## [17] "isDecoy"                   "post"                     
## [19] "pre"                       "end"                      
## [21] "start"                     "accession"                
## [23] "length"                    "description"              
## [25] "pepSeq"                    "modified"                 
## [27] "modification"              "idFile"                   
## [29] "spectrumFile"              "databaseFile"             
## [31] "peptide"

Here, we summarise a few steps and redirect the reader to the package’s vignette for more details:

9.0.1 Analysis of peptide sequences

Cleaning irregular cleavages at the termini of the peptides and missing cleavage site within the peptide sequences. The following two function call create the new numMisCleavages and numIrrCleabages columns in the MSnID object

msnid <- assess_termini(msnid, validCleavagePattern="[KR]\\.[^P]")
msnid <- assess_missed_cleavages(msnid, missedCleavagePattern="[KR](?=[^P$])")

9.1 Trimming the data

Now, we can use the apply_filter function to effectively apply filters. The strings passed to the function represent expressions that will be evaludated, this keeping only PSMs that have 0 irregular cleavages and 2 or less missed cleavages.

msnid <- apply_filter(msnid, "numIrregCleavages == 0")
msnid <- apply_filter(msnid, "numMissCleavages <= 2")
show(msnid)
## MSnID object
## Working directory: "."
## #Spectrum Files:  1 
## #PSMs: 7838 at 17 % FDR
## #peptides: 5598 at 23 % FDR
## #accessions: 3759 at 53 % FDR

9.2 Parent ion mass errors

Using "calculatedMassToCharge" and "experimentalMassToCharge", the mass_measurement_error function calculates the parent ion mass measurement error in parts per million.

summary(mass_measurement_error(msnid))
##       Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
## -2184.0640    -0.6992     0.0000    17.6146     0.7512  2012.5178

We then filter any matches that do not fit the +/- 20 ppm tolerance

msnid <- apply_filter(msnid, "abs(mass_measurement_error(msnid)) < 20")
summary(mass_measurement_error(msnid))
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -19.7797  -0.5866   0.0000  -0.2970   0.5713  19.6758

9.3 Filtering criteria

Filtering of the identification data will rely on

  • -log10 transformed MS-GF+ Spectrum E-value, reflecting the goodness of match experimental and theoretical fragmentation patterns
msnid$msmsScore <- -log10(msnid$`MS-GF:SpecEValue`)
  • the absolute mass measurement error (in ppm units) of the parent ion
msnid$absParentMassErrorPPM <- abs(mass_measurement_error(msnid))

MS2 fiters are handled by a special MSnIDFilter class objects, where individual filters are set by name (that is present in names(msnid)) and comparison operator (>, <, = , …) defining if we should retain hits with higher or lower given the threshold and finally the threshold value itself.

filtObj <- MSnIDFilter(msnid)
filtObj$absParentMassErrorPPM <- list(comparison="<", threshold=10.0)
filtObj$msmsScore <- list(comparison=">", threshold=10.0)
show(filtObj)
## MSnIDFilter object
## (absParentMassErrorPPM < 10) & (msmsScore > 10)

We can then evaluate the filter on the identification data object, which return the false discovery rate and number of retained identifications for the filtering criteria at hand.

evaluate_filter(msnid, filtObj)
##           fdr    n
## PSM         0 3807
## peptide     0 2455
## accession   0 1009

9.4 Filter optimisation

Rather than setting filtering values by hand, as shown above, these can be set automativally to meet a specific false discovery rate.

filtObj.grid <- optimize_filter(filtObj, msnid, fdr.max=0.01,
                                method="Grid", level="peptide",
                                n.iter=500)
show(filtObj.grid)
## MSnIDFilter object
## (absParentMassErrorPPM < 3) & (msmsScore > 7.4)
evaluate_filter(msnid, filtObj.grid)
##                   fdr    n
## PSM       0.004097561 5146
## peptide   0.006447651 3278
## accession 0.021996616 1208

Filters can eventually be applied (rather than just evaluated) using the apply_filter function.

msnid <- apply_filter(msnid, filtObj.grid)
show(msnid)
## MSnID object
## Working directory: "."
## #Spectrum Files:  1 
## #PSMs: 5146 at 0.41 % FDR
## #peptides: 3278 at 0.64 % FDR
## #accessions: 1208 at 2.2 % FDR

And finally, identifications that matched decoy and contaminant protein sequences are removed

msnid <- apply_filter(msnid, "isDecoy == FALSE") 
msnid <- apply_filter(msnid, "!grepl('Contaminant',accession)")
show(msnid)
## MSnID object
## Working directory: "."
## #Spectrum Files:  1 
## #PSMs: 5117 at 0 % FDR
## #peptides: 3251 at 0 % FDR
## #accessions: 1179 at 0 % FDR

The resulting filtered identification data can be exported to a data.frame or to a dedicated MSnSet data structure for quantitative MS data, described below, and further processed and analyses using appropriate statistical tests.

10 High-level data interface

The above sections introduced low-level interfaces to raw and identification results. The MSnbase package provides abstractions for raw data through the MSnExp class and containers for quantification data via the MSnSet class. Both store

  1. the actual assay data (spectra or quantitation matrix, see below), accessed with spectra (or the [, [[ operators) or exprs;
  2. sample metadata, accessed as a data.frame with pData;
  3. feature metadata, accessed as a data.frame with fData.

Another useful slot is processingData, accessed with processingData(.), that records all the processing that objects have undergone since their creation (see examples below).

The readMSData will parse the raw data, extract the MS2 spectra (by default) and construct an MS experiment object of class MSnExp.

(Note that while readMSData supports MS1 data, this is currently not convenient as all the data is read into memory.)

library("MSnbase")
rawFile <- dir(system.file(package = "MSnbase", dir = "extdata"),
               full.name = TRUE, pattern = "mzXML$")
basename(rawFile)
## [1] "dummyiTRAQ.mzXML"
msexp <- readMSData(rawFile, verbose = FALSE, centroided = FALSE)
msexp
## MSn experiment data ("MSnExp")
## Object size in memory: 0.18 Mb
## - - - Spectra data - - -
##  MS level(s): 2 
##  Number of spectra: 5 
##  MSn retention times: 25:1 - 25:2 minutes
## - - - Processing information - - -
## Data loaded: Sun Jul 23 07:53:24 2017 
##  MSnbase version: 2.2.0 
## - - - Meta data  - - -
## phenoData
##   rowNames: dummyiTRAQ.mzXML
##   varLabels: sampleNames
##   varMetadata: labelDescription
## Loaded from:
##   dummyiTRAQ.mzXML 
## protocolData: none
## featureData
##   featureNames: X1.1 X2.1 ... X5.1 (5 total)
##   fvarLabels: spectrum
##   fvarMetadata: labelDescription
## experimentData: use 'experimentData(object)'

MS2 spectra can be extracted as a list of Spectrum2 objects with the spectra accessor or as a subset of the original MSnExp data with the [ operator. Individual spectra can be accessed with [[.

length(msexp)
## [1] 5
msexp[1:2]
## MSn experiment data ("MSnExp")
## Object size in memory: 0.07 Mb
## - - - Spectra data - - -
##  MS level(s): 2 
##  Number of spectra: 2 
##  MSn retention times: 25:1 - 25:2 minutes
## - - - Processing information - - -
## Data loaded: Sun Jul 23 07:53:24 2017 
## Data [numerically] subsetted 2 spectra: Sun Jul 23 07:53:24 2017 
##  MSnbase version: 2.2.0 
## - - - Meta data  - - -
## phenoData
##   rowNames: dummyiTRAQ.mzXML
##   varLabels: sampleNames
##   varMetadata: labelDescription
## Loaded from:
##   dummyiTRAQ.mzXML 
## protocolData: none
## featureData
##   featureNames: X1.1 X2.1
##   fvarLabels: spectrum
##   fvarMetadata: labelDescription
## experimentData: use 'experimentData(object)'
msexp[[2]]
## Object of class "Spectrum2"
##  Precursor: 546.9586 
##  Retention time: 25:2 
##  Charge: 3 
##  MSn level: 2 
##  Peaks count: 1012 
##  Total ion count: 56758067

The identification results stemming from the same raw data file can then be used to add PSM matches.

fData(msexp)
##      spectrum
## X1.1        1
## X2.1        2
## X3.1        3
## X4.1        4
## X5.1        5
## find path to a mzIdentML file
identFile <- dir(system.file(package = "MSnbase", dir = "extdata"),
                 full.name = TRUE, pattern = "dummyiTRAQ.mzid")
basename(identFile)
## [1] "dummyiTRAQ.mzid"
msexp <- addIdentificationData(msexp, identFile)
fData(msexp)
##      spectrum scan number(s) passthreshold rank calculatedmasstocharge
## X1.1        1              1          TRUE    1               645.0375
## X2.1        2              2          TRUE    1               546.9633
## X3.1        3             NA            NA   NA                     NA
## X4.1        4             NA            NA   NA                     NA
## X5.1        5              5          TRUE    1               437.2997
##      experimentalmasstocharge chargestate ms-gf:denovoscore ms-gf:evalue
## X1.1                 645.3741           3                77     79.36958
## X2.1                 546.9586           3                39     13.46615
## X3.1                       NA          NA                NA           NA
## X4.1                       NA          NA                NA           NA
## X5.1                 437.8040           2                 5    366.38422
##      ms-gf:rawscore ms-gf:specevalue assumeddissociationmethod
## X1.1            -39     5.527468e-05                       CID
## X2.1            -30     9.399048e-06                       CID
## X3.1             NA               NA                      <NA>
## X4.1             NA               NA                      <NA>
## X5.1            -42     2.577830e-04                       CID
##      isotopeerror isdecoy post  pre end start       accession length
## X1.1            1   FALSE    A    R 186   170 ECA0984;ECA3829    231
## X2.1            0   FALSE    A    K  62    50         ECA1028    275
## X3.1         <NA>      NA <NA> <NA>  NA    NA            <NA>     NA
## X4.1         <NA>      NA <NA> <NA>  NA    NA            <NA>     NA
## X5.1            1   FALSE    L    K  28    22         ECA0510    166
##                                                                      description
## X1.1 DNA mismatch repair protein;acetolactate synthase isozyme III large subunit
## X2.1          2,3,4,5-tetrahydropyridine-2,6-dicarboxylate N-succinyltransferase
## X3.1                                                                        <NA>
## X4.1                                                                        <NA>
## X5.1                    putative capsular polysacharide biosynthesis transferase
##                 pepseq modified modification          idFile
## X1.1 VESITARHGEVLQLRPK    FALSE           NA dummyiTRAQ.mzid
## X2.1     IDGQWVTHQWLKK    FALSE           NA dummyiTRAQ.mzid
## X3.1              <NA>       NA           NA            <NA>
## X4.1              <NA>       NA           NA            <NA>
## X5.1           LVILLFR    FALSE           NA dummyiTRAQ.mzid
##                  databaseFile nprot npep.prot npsm.prot npsm.pep
## X1.1 erwinia_carotovora.fasta     2         1         1        1
## X2.1 erwinia_carotovora.fasta     1         1         1        1
## X3.1                     <NA>    NA        NA        NA       NA
## X4.1                     <NA>    NA        NA        NA       NA
## X5.1 erwinia_carotovora.fasta     1         1         1        1

The readMSData and addIdentificationData make use of mzR and mzID packages to access the raw and identification data.

Spectra and (parts of) experiments can be extracted and plotted.

msexp[[1]]
## Object of class "Spectrum2"
##  Precursor: 645.3741 
##  Retention time: 25:1 
##  Charge: 3 
##  MSn level: 2 
##  Peaks count: 2921 
##  Total ion count: 668170086
plot(msexp[[1]], full=TRUE)

msexp[1:3]
## MSn experiment data ("MSnExp")
## Object size in memory: 0.1 Mb
## - - - Spectra data - - -
##  MS level(s): 2 
##  Number of spectra: 3 
##  MSn retention times: 25:1 - 25:2 minutes
## - - - Processing information - - -
## Data loaded: Sun Jul 23 07:53:24 2017 
## Data [numerically] subsetted 3 spectra: Sun Jul 23 07:53:26 2017 
##  MSnbase version: 2.2.0 
## - - - Meta data  - - -
## phenoData
##   rowNames: dummyiTRAQ.mzXML
##   varLabels: sampleNames
##   varMetadata: labelDescription
## Loaded from:
##   dummyiTRAQ.mzXML 
## protocolData: none
## featureData
##   featureNames: X1.1 X2.1 X3.1
##   fvarLabels: spectrum scan number(s) ... npsm.pep (30 total)
##   fvarMetadata: labelDescription
## experimentData: use 'experimentData(object)'
plot(msexp[1:3], full=TRUE)

11 Quantitative proteomics

There are a wide range of proteomics quantitation techniques that can broadly be classified as labelled vs. label-free, depending whether the features are labelled prior the MS acquisition and the MS level at which quantitation is inferred, namely MS1 or MS2.

Label-free Labelled
MS1 XIC SILAC, 15N
MS2 Counting iTRAQ, TMT

In terms of raw data quantitation, most efforts have been devoted to MS2-level quantitation. Label-free XIC quantitation has however been addressed in the frame of metabolomics data processing by the xcms infrastructure.

An MSnExp is converted to an MSnSet by the quantitation method. Below, we use the iTRAQ 4-plex isobaric tagging strategy (defined by the iTRAQ4 parameter; other tags are available) and the trapezoidation method to calculate the area under the isobaric reporter peaks.

plot(msexp[[1]], full=TRUE, reporters = iTRAQ4)

msset <- quantify(msexp, method = "trap", reporters = iTRAQ4, verbose=FALSE)
exprs(msset)
##      iTRAQ4.114 iTRAQ4.115 iTRAQ4.116 iTRAQ4.117
## X1.1   4483.320   4873.996   6743.441   4601.378
## X2.1   1918.082   1418.040   1117.601   1581.954
## X3.1  15210.979  15296.256  15592.760  16550.502
## X4.1   4133.103   5069.983   4724.845   4694.801
## X5.1  11947.881  13061.875  12809.491  12911.479
processingData(msset)
## - - - Processing information - - -
## Data loaded: Sun Jul 23 07:53:24 2017 
## iTRAQ4 quantification by trapezoidation: Sun Jul 23 07:53:30 2017 
##  MSnbase version: 2.2.0

Other MS2 quantitation methods available in quantify include the (normalised) spectral index SI and (normalised) spectral abundance factor SAF or simply a simple count method.

exprs(si <- quantify(msexp, method = "SIn"))     
##         dummyiTRAQ.mzXML
## ECA0510      0.003588641
## ECA1028      0.001470129
exprs(saf <- quantify(msexp, method = "NSAF"))
##         dummyiTRAQ.mzXML
## ECA0510        0.6235828
## ECA1028        0.3764172

Note that spectra that have not been assigned any peptide (NA) or that match non-unique peptides (npsm > 1) are discarded in the counting process.

See also The isobar package supports quantitation from centroided mgf peak lists or its own tab-separated files that can be generated from Mascot and Phenyx vendor files.

12 Importing third-party quantitation data

The PSI mzTab file format is aimed at providing a simpler (than XML formats) and more accessible file format to the wider community. It is composed of a key-value metadata section and peptide/protein/small molecule tabular sections.

Note that below, we specify version 0.9 (that generates the warning) to fit with the file. For recent files, the version argument would be ignored to use the recent importer.

mztf <- pxget(px, "F063721.dat-mztab.txt")
## Downloading 1 file
(mzt <- readMzTabData(mztf, what = "PEP", version = "0.9"))
## Warning: Version 0.9 is deprecated. Please see '?readMzTabData' and '?
## MzTab' for details.
## MSnSet (storageMode: lockedEnvironment)
## assayData: 1528 features, 6 samples 
##   element names: exprs 
## protocolData: none
## phenoData
##   sampleNames: sub[1] sub[2] ... sub[6] (6 total)
##   varLabels: abundance
##   varMetadata: labelDescription
## featureData
##   featureNames: 1 2 ... 1528 (1528 total)
##   fvarLabels: sequence accession ... uri (14 total)
##   fvarMetadata: labelDescription
## experimentData: use 'experimentData(object)'
## Annotation:  
## - - - Processing information - - -
## mzTab read: Sun Jul 23 07:53:35 2017 
##  MSnbase version: 2.2.0

It is also possible to import arbitrary spreadsheets as MSnSet objects into R with the readMSnSet2 function. The main 2 arguments of the function are (1) a text-based spreadsheet and (2) column names of indices that identify the quantitation data. The latter can be queried with the getEcols function.

csv <- dir(system.file ("extdata" , package = "pRolocdata"),
           full.names = TRUE, pattern = "pr800866n_si_004-rep1.csv")
getEcols(csv, split = ",")
##  [1] "\"Protein ID\""              "\"FBgn\""                   
##  [3] "\"Flybase Symbol\""          "\"No. peptide IDs\""        
##  [5] "\"Mascot score\""            "\"No. peptides quantified\""
##  [7] "\"area 114\""                "\"area 115\""               
##  [9] "\"area 116\""                "\"area 117\""               
## [11] "\"PLS-DA classification\""   "\"Peptide sequence\""       
## [13] "\"Precursor ion mass\""      "\"Precursor ion charge\""   
## [15] "\"pd.2013\""                 "\"pd.markers\""
ecols <- 7:10
res <- readMSnSet2(csv, ecols)
head(exprs(res))
##   area.114 area.115 area.116 area.117
## 1 0.379000 0.281000 0.225000 0.114000
## 2 0.420000 0.209667 0.206111 0.163889
## 3 0.187333 0.167333 0.169667 0.476000
## 4 0.247500 0.253000 0.320000 0.179000
## 5 0.216000 0.183000 0.342000 0.259000
## 6 0.072000 0.212333 0.573000 0.142667
head(fData(res))
##   Protein.ID        FBgn Flybase.Symbol No..peptide.IDs Mascot.score
## 1    CG10060 FBgn0001104    G-ialpha65A               3       179.86
## 2    CG10067 FBgn0000044         Act57B               5       222.40
## 3    CG10077 FBgn0035720        CG10077               5       219.65
## 4    CG10079 FBgn0003731           Egfr               2        86.39
## 5    CG10106 FBgn0029506        Tsp42Ee               1        52.10
## 6    CG10130 FBgn0010638      Sec61beta               2        79.90
##   No..peptides.quantified PLS.DA.classification Peptide.sequence
## 1                       1                    PM                 
## 2                       9                    PM                 
## 3                       3                                       
## 4                       2                    PM                 
## 5                       1                              GGVFDTIQK
## 6                       3              ER/Golgi                 
##   Precursor.ion.mass Precursor.ion.charge     pd.2013 pd.markers
## 1                                                  PM    unknown
## 2                                                  PM    unknown
## 3                                             unknown    unknown
## 4                                                  PM    unknown
## 5            626.887                    2 Phenotype 1    unknown
## 6                                            ER/Golgi         ER

13 Data processing and analysis

13.1 Raw data processing

For raw data processing look at MSnbase’s clean, smooth, pickPeaks, removePeaks and trimMz for MSnExp and spectra processing methods.

The MALDIquantand xcms packages also features a wide range of raw data processing methods on their own ad hoc data instance types.

13.2 Processing and normalisation

Each different types of quantitative data will require their own pre-processing and normalisation steps. Both isobar and MSnbase allow to correct for isobaric tag impurities normalise the quantitative data.

data(itraqdata)
qnt <- quantify(itraqdata, method = "trap",
                reporters = iTRAQ4, verbose = FALSE)
impurities <- matrix(c(0.929,0.059,0.002,0.000,
                       0.020,0.923,0.056,0.001,
                       0.000,0.030,0.924,0.045,
                       0.000,0.001,0.040,0.923),
                     nrow=4, byrow = TRUE)
## or, using makeImpuritiesMatrix()
## impurities <- makeImpuritiesMatrix(4)
qnt.crct <- purityCorrect(qnt, impurities)
processingData(qnt.crct)
## - - - Processing information - - -
## Data loaded: Wed May 11 18:54:39 2011 
## Updated from version 0.3.0 to 0.3.1 [Fri Jul  8 20:23:25 2016] 
## iTRAQ4 quantification by trapezoidation: Sun Jul 23 07:53:39 2017 
## Purity corrected: Sun Jul 23 07:53:39 2017 
##  MSnbase version: 1.1.22

Various normalisation methods can be applied the MSnSet instances using the normalise method: variance stabilisation (vsn), quantile (quantiles), median or mean centring (center.media or center.mean), …

qnt.crct.nrm <- normalise(qnt.crct, "quantiles") 

The combineFeatures method combines spectra/peptides quantitation values into protein data. The grouping is defined by the groupBy parameter, which is generally taken from the feature metadata (protein accessions, for example).

## arbitraty grouping
g <- factor(c(rep(1, 25), rep(2, 15), rep(3, 15)))
g
##  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2
## [36] 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
## Levels: 1 2 3
prt <- combineFeatures(qnt.crct.nrm, groupBy = g, fun = "sum")
processingData(prt)
## - - - Processing information - - -
## Data loaded: Wed May 11 18:54:39 2011 
## Updated from version 0.3.0 to 0.3.1 [Fri Jul  8 20:23:25 2016] 
## iTRAQ4 quantification by trapezoidation: Sun Jul 23 07:53:39 2017 
## Purity corrected: Sun Jul 23 07:53:39 2017 
## Normalised (quantiles): Sun Jul 23 07:53:39 2017 
## Combined 55 features into 3 using sum: Sun Jul 23 07:53:39 2017 
##  MSnbase version: 2.2.0

Finally, proteomics data analysis is generally hampered by missing values. Missing data imputation is a sensitive operation whose success will be guided by many factors, such as degree and (non-)random nature of the missingness.

Below, missing values are randomly assigned to our test data and visualised on a heatmap.

set.seed(1)
qnt0 <- qnt
exprs(qnt0)[sample(prod(dim(qnt0)), 10)] <- NA
table(is.na(qnt0))
## 
## FALSE  TRUE 
##   209    11
image(qnt0)

Missing value in MSnSet instances can be filtered out and imputed using the filterNA and impute functions.

## remove features with missing values
qnt00 <- filterNA(qnt0)
dim(qnt00)
## [1] 44  4
any(is.na(qnt00))
## [1] FALSE
## impute missing values using knn imputation
qnt.imp <- impute(qnt0, method = "knn")
dim(qnt.imp)
## [1] 55  4
any(is.na(qnt.imp))
## [1] FALSE

There are various methods to perform data imputation, as described in ?impute.

14 Statistical analysis

R in general and Bioconductor in particular are well suited for the statistical analysis of data. Several packages provide dedicated resources for proteomics data:

  • MSstats: A set of tools for statistical relative protein significance analysis in DDA, SRM and DIA experiments. Data stored in data.frame or MSnSet objects can be used as input.

  • msmsTests: Statistical tests for label-free LC-MS/MS data by spectral counts, to discover differentially expressed proteins between two biological conditions. Three tests are available: Poisson GLM regression, quasi-likelihood GLM regression, and the negative binomial of the edgeR package. All can be readily applied on MSnSet instances produced, for example by MSnID.

  • isobar also provides dedicated infrastructure for the statistical analysis of isobaric data.

15 Machine learning

The MLInterfaces package provides a unified interface to a wide range of machine learning algorithms. Initially developed for microarray and ExpressionSet instances, the pRoloc package enables application of these algorithms to MSnSet data.

15.1 Classification

The example below uses knn with the 5 closest neighbours as an illustration to classify proteins of unknown sub-cellular localisation to one of 9 possible organelles.

library("MLInterfaces")
library("pRoloc")
library("pRolocdata")
data(dunkley2006)
traininds <- which(fData(dunkley2006)$markers != "unknown")
ans <- MLearn(markers ~ ., data = t(dunkley2006), knnI(k = 5), traininds)
ans
## MLInterfaces classification output container
## The call was:
## MLearn(formula = markers ~ ., data = t(dunkley2006), .method = knnI(k = 5), 
##     trainInd = traininds)
## Predicted outcome distribution for test set:
## 
##      ER lumen   ER membrane         Golgi Mitochondrion       Plastid 
##             5           140            67            51            29 
##            PM      Ribosome           TGN       vacuole 
##            89            31             6            10 
## Summary of scores on test set (use testScores() method for details):
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.4000  1.0000  1.0000  0.9332  1.0000  1.0000

15.2 Clustering

15.2.1 kmeans

kcl <- MLearn( ~ ., data = dunkley2006, kmeansI, centers = 12)
kcl
## clusteringOutput: partition table
## 
##   1   2   3   4   5   6   7   8   9  10  11  12 
##  55  28  29  18  60  30  44 145  78 116  33  53 
## The call that created this object was:
## MLearn(formula = ~., data = dunkley2006, .method = kmeansI, centers = 12)
plot(kcl, exprs(dunkley2006))

A wide range of classification and clustering algorithms are also available, as described in the ?MLearn documentation page. The pRoloc package also uses MSnSet instances as input and ,while being conceived with the analysis of spatial/organelle proteomics data in mind, is applicable many use cases.

16 Annotation

All the Bioconductor annotation infrastructure, such as biomaRt, GO.db, organism specific annotations, .. are directly relevant to the analysis of proteomics data. A total of 193 ontologies, including some proteomics-centred annotations such as the PSI Mass Spectrometry Ontology, Molecular Interaction (PSI MI 2.5) or Protein Modifications are available through the rols

library("rols")
res <- OlsSearch(q = "ESI", ontology = "MS", exact = TRUE)
res
## Object of class 'OlsSearch':
##   ontolgy: MS 
##   query: ESI 
##   requested: 20 (out of 1)
##   response(s): 0

There is a single exact match (default is to retrieve 20 results), that can be retrieved and coreced to a Terms or data.frame object with

res <- olsSearch(res)
as(res, "Terms")
## Object of class 'Terms' with 1 entries
##  From the MS ontology
## MS:1000073
as(res, "data.frame")
##                                                   id
## 1 ms:class:http://purl.obolibrary.org/obo/MS_1000073
##                                         iri short_form     obo_id
## 1 http://purl.obolibrary.org/obo/MS_1000073 MS_1000073 MS:1000073
##                     label
## 1 electrospray ionization
##                                                                                                                                                                                                                                                                                                                                                                                                                                                  description
## 1 A process in which ionized species in the gas phase are produced from an analyte-containing solution via highly charged fine droplets, by means of spraying the solution from a narrow-bore needle tip at atmospheric pressure in the presence of a high electric field. When a pressurized gas is used to aid in the formation of a stable spray, the term pneumatically assisted electrospray ionization is used. The term ion spray is not recommended.
##   ontology_name ontology_prefix  type is_defining_ontology
## 1            ms              MS class                 TRUE

Data from the Human Protein Atlas is available via the hpar package.

17 Other relevant packages/pipelines

  • Analysis of post translational modification with isobar.
  • Analysis of label-free data from a Synapt G2 (including ion mobility) with synapter.
  • Analysis of spatial proteomics data with pRoloc.
  • Analysis of MALDI data with the MALDIquant package.
  • Access to the Proteomics Standard Initiative Common QUery InterfaCe with the PSICQUIC package.

Additional relevant packages are described in the RforProteomics vignettes.

18 Session information

## R version 3.4.1 (2017-06-30)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Debian GNU/Linux 9 (stretch)
## 
## Matrix products: default
## BLAS/LAPACK: /usr/lib/libopenblasp-r0.2.19.so
## 
## attached base packages:
## [1] stats4    parallel  methods   stats     graphics  grDevices utils    
## [8] datasets  base     
## 
## other attached packages:
##  [1] hpar_1.18.1           rols_2.4.0            MSGFplus_1.10.1      
##  [4] pRolocdata_1.14.0     pRoloc_1.16.1         MLInterfaces_1.56.0  
##  [7] cluster_2.0.6         annotate_1.54.0       XML_3.98-1.9         
## [10] AnnotationDbi_1.38.1  IRanges_2.10.2        S4Vectors_0.14.3     
## [13] rpx_1.12.1            MSnID_1.10.0          mzID_1.14.0          
## [16] nloptr_1.0.4          RforProteomics_1.14.0 MSnbase_2.2.0        
## [19] ProtGenerics_1.8.0    BiocParallel_1.10.1   mzR_2.10.0           
## [22] Rcpp_0.12.11          Biobase_2.36.2        BiocGenerics_0.22.0  
## [25] BiocStyle_2.4.1       BiocInstaller_1.26.0  shiny_1.0.3          
## [28] rmarkdown_1.6         knitr_1.16           
## 
## loaded via a namespace (and not attached):
##   [1] questionr_0.6.1               R.utils_2.5.0                
##   [3] RUnit_0.4.31                  lme4_1.1-13                  
##   [5] RSQLite_2.0                   htmlwidgets_0.9              
##   [7] grid_3.4.1                    trimcluster_0.1-2            
##   [9] lpSolve_5.6.13                rda_1.0.2-2                  
##  [11] munsell_0.4.3                 codetools_0.2-15             
##  [13] preprocessCore_1.38.1         miniUI_0.1.1                 
##  [15] colorspace_1.3-2              Category_2.42.1              
##  [17] highr_0.6                     rstudioapi_0.6               
##  [19] robustbase_0.92-7             labeling_0.3                 
##  [21] hwriter_1.3.2                 bit64_0.9-7                  
##  [23] ggvis_0.4.3                   rprojroot_1.2                
##  [25] randomForest_4.6-12           diptest_0.75-7               
##  [27] R6_2.2.2                      doParallel_1.0.10            
##  [29] gridSVG_1.5-1                 rmdformats_0.3.3             
##  [31] flexmix_2.3-14                bitops_1.0-6                 
##  [33] assertthat_0.2.0              scales_0.4.1                 
##  [35] nnet_7.3-12                   gtable_0.2.0                 
##  [37] affy_1.54.0                   biocViews_1.44.0             
##  [39] rlang_0.1.1                   MatrixModels_0.4-1           
##  [41] genefilter_1.58.1             splines_3.4.1                
##  [43] lazyeval_0.2.0                ModelMetrics_1.1.0           
##  [45] impute_1.50.1                 hexbin_1.27.1                
##  [47] yaml_2.1.14                   reshape2_1.4.2               
##  [49] threejs_0.2.2                 backports_1.1.0              
##  [51] httpuv_1.3.5                  RBGL_1.52.0                  
##  [53] caret_6.0-76                  tools_3.4.1                  
##  [55] bookdown_0.4                  ggplot2_2.2.1                
##  [57] affyio_1.46.0                 RColorBrewer_1.1-2           
##  [59] proxy_0.4-17                  plyr_1.8.4                   
##  [61] base64enc_0.1-3               progress_1.1.2               
##  [63] zlibbioc_1.22.0               RCurl_1.95-4.8               
##  [65] prettyunits_1.0.2             rpart_4.1-11                 
##  [67] viridis_0.4.0                 sampling_2.8                 
##  [69] sfsmisc_1.1-1                 magrittr_1.5                 
##  [71] data.table_1.10.4             SparseM_1.77                 
##  [73] pcaMethods_1.68.0             mvtnorm_1.0-6                
##  [75] whisker_0.3-2                 R.cache_0.12.0               
##  [77] mime_0.5                      evaluate_0.10.1              
##  [79] xtable_1.8-2                  pbkrtest_0.4-7               
##  [81] mclust_5.3                    gridExtra_2.2.1              
##  [83] compiler_3.4.1                biomaRt_2.32.1               
##  [85] tibble_1.3.3                  minqa_1.2.4                  
##  [87] R.oo_1.21.0                   htmltools_0.3.6              
##  [89] mgcv_1.8-17                   DBI_0.7                      
##  [91] MASS_7.3-47                   fpc_2.1-10                   
##  [93] Matrix_1.2-10                 car_2.1-5                    
##  [95] vsn_3.44.0                    R.methodsS3_1.7.1            
##  [97] gdata_2.18.0                  mlbench_2.1-1                
##  [99] bindr_0.1                     pkgconfig_2.0.1              
## [101] MALDIquant_1.16.2             xml2_1.1.1                   
## [103] foreach_1.4.3                 stringr_1.2.0                
## [105] digest_0.6.12                 pls_2.6-0                    
## [107] graph_1.54.0                  dendextend_1.5.2             
## [109] GSEABase_1.38.0               curl_2.7                     
## [111] kernlab_0.9-25                gtools_3.5.0                 
## [113] quantreg_5.33                 modeltools_0.2-21            
## [115] nlme_3.1-131                  jsonlite_1.5                 
## [117] bindrcpp_0.2                  interactiveDisplay_1.14.0    
## [119] viridisLite_0.2.0             limma_3.32.3                 
## [121] lattice_0.20-35               httr_1.2.1                   
## [123] DEoptimR_1.0-8                survival_2.41-3              
## [125] interactiveDisplayBase_1.14.0 glue_1.1.1                   
## [127] FNN_1.1                       gbm_2.1.3                    
## [129] prabclus_2.2-6                iterators_1.0.8              
## [131] bit_1.1-12                    class_7.3-14                 
## [133] stringi_1.1.5                 blob_1.1.0                   
## [135] memoise_1.1.0                 dplyr_0.7.1                  
## [137] e1071_1.6-8