This material available under a creative common
CC-BY license. You
are free to share (copy and redistribute the material in any
medium or format) and adapt (remix, transform, and build upon
the material) for any purpose, even commercially.
Biological information flow
Information flow in biological systems (Source Central dogma of biology on Wikipedia).
There are three main biological entities that respectively store
information, act as data intermediates, and the functional units, and
the figure below show how information flows between these three
levels.
DNA, that lives in the nucleus of cells, is the central
information storage mechanism, and encodes the blueprint of the
functional units as genes. DNA is transcribed into messenger
RNA (mRNA), that re-localises outside the nucleus and is further
processed into its mature exon-only form after removal of the
non-coding introns sequences. Finally, the mRNA is translated by the
ribosomial machinery into proteins directly into the endoplasmic
reticulum (ER) where they are then redirected to their final
destination.
Sub-cellular structure of an animal cell (Source Cell biology on Wikipedia).
In addition to the standard information worflow where DNA is
transcribed into RNA that itself is translated into proteins,
there is also reverse transcription, that generates
(complementary) DNA from and RNA molecule, as well as replication of
DNA (during cell division) and RNA molecules.
Why studying proteins
Proteins as the functional units in all living organisms, and they are
highly dynamic. The caterpillar and the resulting butterfly have the
same genome. The complement of all the expressed proteins, termed the
proteome is however very different.
The metamorphosis from a caterpilar to a monarch butterfly. (Image from Phys.org)
There are different modalities of the proteome that are of
interest. In addition to the presence or amount of protein in a
biological samples, it is also important to study the interactions
between proteins forming protein-protein complexes, the presence of
post-transcriptional modification (such as, for example,
phosphorylations), the rate at which proteins are produced and
degraded, or where the proteins reside inside a cell.
The technique of choice to study proteins in a high throughput way is
mass spectrometry.
A beginner’s guide to mass spectrometry–based proteomics(Sinha and Mann 2020) is an approachable introduction to sample preparation,
mass spectrometry and data analysis.
Setup
We are going to use Bioconductor (Huber et al. 2015) packages from the R
for Mass Spectrometry
initiative. The associated
book provide further
details and state-of-the-art infrastructure for MS and proteomics data
analysis.
The R for Mass Spectrometry intiative.
Packages can be install with the BiocManager package, available from
CRAN. If BiocManager isn’t available on your computer, install it
with:
install.packages("BiocManager")
Any packages (such as Spectra) and their dependencies can be
installed with:
BiocManager::install("Spectra")
For additional information on how to analyse mass spectrometry-based
proteomics data, refer to (L. Gatto and Christoforou 2014) and (Laurent Gatto 2019), or explore
the
proteomics-
and mass
spectrometry-related
packages on the Bioconductor page
How does mass spectrometry work?
Mass spectrometry (MS) is a technology that separates charged
molecules (ions) based on their mass to charge ratio (M/Z). It is
often coupled to chromatography (liquid LC, but can also be gas-based
GC). The time an analytes takes to elute from the chromatography
column is the retention time.
A chromatogram, illustrating the total amount of analytes over the retention time.
An mass spectrometer is composed of three components:
The source, that ionises the molecules: examples are Matrix-assisted
laser desorption/ionisation (MALDI) or electrospray ionisation.
(ESI)
The analyser, that separates the ions: Time of flight (TOF) or Orbitrap.
The detector that quantifies the ions.
When using mass spectrometry for proteomics, the proteins are first
digested with a protease such as trypsin. In mass shotgun proteomics,
the analytes assayed in the mass spectrometer are peptides.
Often, ions are subjected to more than a single MS round. After a
first round of separation, the peaks in the spectra, called MS1
spectra, represent peptides. At this stage, the only information we
possess about these peptides are their retention time and their
mass-to-charge (we can also infer their charge be inspecting their
isotopic envelope, i.e the peaks of the individual isotopes, see
below), which is not enough to infer their identify (i.e. their
sequence).
In MSMS (or MS2), the settings of the mass spectrometer are set
automatically to select a certain number of MS1 peaks (for example
20). Once a narrow M/Z range has been selected (corresponding to one
high-intensity peak, a peptide, and some background noise), it is
fragmented (using for example collision-induced dissociation (CID),
higher energy collisional dissociation (HCD) or electron-transfer
dissociation (ETD)). The fragment ions are then themselves separated
in the analyser to produce a MS2 spectrum. The unique fragment ion
pattern can then be used to infer the peptide sequence using de novo
sequencing (when the spectrum is of high enough quality) of using a
search engine such as, for example Mascot, MSGF+, …, that will match
the observed, experimental spectrum to theoratical spectra (see
details below).
Schematics of a mass spectrometer and two rounds of MS.
The animation below show how 25 ions different ions (i.e. having
different M/Z values) are separated throughout the MS analysis and are
eventually detected (i.e. quantified). The final frame shows the
hypothetical spectrum.
Separation and detection of ions in a mass spectrometer.
The figures below illustrate the two rounds of MS. The spectrum on the
left is an MS1 spectrum acquired after 21 minutes and 3 seconds of
elution. 10 peaks, highlited by dotted vertical lines, were selected
for MS2 analysis. The peak at M/Z 460.79 (488.8) is highlighted by a
red (orange) vertical line on the MS1 spectrum and the fragment
spectra are shown on the MS2 spectrum on the top (bottom) right
figure.
Parent ions in the MS1 spectrum (left) and two sected fragment ions MS2 spectra (right).
The figures below represent the 3 dimensions of MS data: a set of
spectra (M/Z and intensity) of retention time, as well as the
interleaved nature of MS1 and MS2 (and there could be more levels)
data.
MS1 spectra (blue) over retention time (left). MS2 spectra (pink) interleaved between two MS1 spectra (right),
Reading and accessing MS data
Let’s read a very small raw MS data file into R using the Spectra()
from the Spectra package
(Rainer et al. 2026). The file that we are going to load is available in
the MsDataHub data package.
Load the Spectra package
library("Spectra")
Get the path to the raw MS data file. When you run this command for
the first time, it will download the MS data file in a local cache
so that next time you need it, it can be retrieved directly.
The object that is returned by Spectra() is of class Specta, that
can store, access and manipulate raw MS data. Note that here we are
focusing on MS-based proteomics data, but this also applied to
MS-based metabolomics data.
We can find out how many spectra are available in that data using
the function length. Full MS acquisitions would contain hundreds
of thousands spectra.
length(sp)
## [1] 7534
We can use various accessor function to get the MS level of these
spectra, their retention time, or the M/Z and intensity of the
precursor peaks of the ion corresponding to the MS2 spectra.
## [1] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [26] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [51] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [76] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [ reached 'max' / getOption("max.print") -- omitted 7434 entries ]
precursorIntensity(sp)
## [1] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [26] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [51] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [76] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [ reached 'max' / getOption("max.print") -- omitted 7434 entries ]
We can also extract subsets of the full experiment or individual
spectra using [ and plot them.
sp[3]
## MSn data (Spectra) with 1 spectra in a MsBackendMzR backend:
## msLevel rtime scanIndex
## <integer> <numeric> <integer>
## 1 1 1.8524 3
## ... 34 more variables/columns.
##
## file(s):
## 4189228548251_7858
plotSpectra(sp[1119])
Visualisation of the 3rd MS spectrum in our small test data set.
Exercises
Spectra objects also contains additional metadata for each spectrum,
that can be accessed, as a DataFrame, with spectraData. What
variables are available?
spectraData(sp)
## DataFrame with 7534 rows and 35 columns
## msLevel rtime acquisitionNum scanIndex dataStorage
## <integer> <numeric> <integer> <integer> <character>
## 1 1 0.4584 1 1 /home/lgatto/.cache/..
## dataOrigin centroided smoothed polarity precScanNum
## <character> <logical> <logical> <integer> <integer>
## 1 /home/lgatto/.cache/.. FALSE NA 1 NA
## precursorMz precursorIntensity precursorCharge collisionEnergy
## <numeric> <numeric> <integer> <numeric>
## 1 NA NA NA NA
## isolationWindowLowerMz isolationWindowTargetMz isolationWindowUpperMz
## <numeric> <numeric> <numeric>
## 1 NA NA NA
## peaksCount totIonCurrent basePeakMZ basePeakIntensity electronBeamEnergy
## <integer> <numeric> <numeric> <numeric> <numeric>
## 1 25800 9187505 445.12 984171 NA
## ionisationEnergy lowMZ highMZ mergedScan mergedResultScanNum
## <numeric> <numeric> <numeric> <integer> <integer>
## 1 0 399.999 2008.45 NA NA
## mergedResultStartScanNum mergedResultEndScanNum injectionTime
## <integer> <integer> <numeric>
## 1 NA NA 50.3826
## filterString spectrumId ionMobilityDriftTime
## <character> <character> <numeric>
## 1 FTMS + p NSI Full ms.. controllerType=0 con.. NA
## scanWindowLowerLimit scanWindowUpperLimit
## <numeric> <numeric>
## 1 400 2000
## [ reached 'max' / getOption("max.print") -- omitted 10 rows ]
The actual MS data, i.e. the M/Z and intensity pairs are accessed with
the peaksData() function, that returns a list of matrices. Extract
the peaks of the spectra 1, 1119, 2235 and 7500.
peaksData(sp[c(1, 1119, 2235, 7500)])
## List of length 4
Identification
The raw data is still a long way of obtaining biologically relevant
proteomics data. The first step to obtain proteomics data is to
identify the peptides that have been acquired in the MS. Peptide
identification work by comparing expected and observed spectra. As
shown below, when a precursor peptide ion is fragmented in a CID cell,
it breaks at specific bonds, producing sets of peaks (a, b, c
and x, y, z) that can be predicted.
Peptide fragmentation.
It is thus possible to calculate the expected set of fagment peaks for
a given peptide, such as SIGFEGDSIGR below. We will make use of the
PSMatch package (Deflandre, Gibb, and Gatto 2026) to handle
identification data.
## mz ion type pos z seq peptide
## 1 88.03931 b1 b 1 1 S SIGFEGDSIGR
## 2 201.12337 b2 b 2 1 SI SIGFEGDSIGR
## 3 258.14483 b3 b 3 1 SIG SIGFEGDSIGR
## 4 405.21324 b4 b 4 1 SIGF SIGFEGDSIGR
## 5 534.25583 b5 b 5 1 SIGFE SIGFEGDSIGR
## 6 591.27729 b6 b 6 1 SIGFEG SIGFEGDSIGR
## 7 706.30423 b7 b 7 1 SIGFEGD SIGFEGDSIGR
## 8 793.33626 b8 b 8 1 SIGFEGDS SIGFEGDSIGR
## 9 906.42032 b9 b 9 1 SIGFEGDSI SIGFEGDSIGR
## 10 963.44178 b10 b 10 1 SIGFEGDSIG SIGFEGDSIGR
## 11 175.11895 y1 y 1 1 R SIGFEGDSIGR
## 12 232.14041 y2 y 2 1 GR SIGFEGDSIGR
## 13 345.22447 y3 y 3 1 IGR SIGFEGDSIGR
## 14 432.25650 y4 y 4 1 SIGR SIGFEGDSIGR
## [ reached 'max' / getOption("max.print") -- omitted 18 rows ]
To illustrate this, we are going to use data from the following
file. Once
downloaded, let’s load it to obtain a Spectra object that contains
identified spectra.
spi <- readRDS("itraqdata2-spectra.rds")
The identification process aims to to compare obseved and expected
peaks. If there is a good match, the MS2 spectrum is assigned the
peptide sequence. Below, we use the plotSpectraPTM() function that
automatically annotates the peaks based on the
s <- "SIGFEGDSIGR"
plotSpectraPTM(spi[14])
Matching observed and expected peaks.
It is also possible to plot 2 spectra (that for example match the same
sequence) to compare them directly
spi[c(25, 28)]$sequence
## [1] "IMIDLDGTENK" "IMIDLDGTENK"
plotSpectraMirror(spi[25], spi[28])
Direct comparison of 2 MS2 spectra.
In a full experiment, all possible peptides from the known (or
relevant) proteome of interest (such as databases that can be
downloaded from the UniProt
site1 The Universal Protein Resource (UniProt) is a freely and
accessible comprehensive resource for protein sequence and
annotation data.) are compared to the
millions of observed spectra.
Statistical challenge The matching between millions
of observed and possible spectra causes real challenges due to the large
search space and the risk of false positives. See (Käll et al. 2008) for further reading.
From the list of identified peptides, it is then necessary to infer the
most provable proteins that were present in the biological sample.
Statistical challenge Protein inference is a
difficult task, as peptides often match multiple proteins (either
different isoforms or proteins stemming from different gene but with
identical domains), which leads to the definition of protein groups,
i.e. sets of proteins that can’t be distinguished with the set of
identified peptides at hand. See (Nesvizhskii and Aebersold 2005) for further reading.
Quantitation
The last step of MS data processing is to quantify peptide abundances
in the biological samples. The table below summarises the different
possibilites depending whether the proteins or peptides are labelled,
and whether the quantitation is performed in MS1 or MS2.
Label-free
Labelled
MS1
XIC
SILAC, 15N
MS2
Counting
iTRAQ, TMT
Label-free MS2: Spectral counting
In spectral counting, on simply counts the number of quantified
peptides that are assigned to a protein.
Labelled MS2: Isobaric tagging
Isobaric tagging refers to the labelling using isobaric tags,
i.e. chemical tags that have the same mass and hence can’t be
distinguish by the spectrometer. The peptides of different samples (4,
6, 10 or 11) are labelled with different tags and combined prior to
mass spectrometry acquisition. Given that they are isobaric, all
identical peptides, irrespective of the tag and this the sample of
origin, are co-analysed, up to fragmentation prior to MS2
analysis. During fragmentation, the isobaric tags fall of, fragment
themselves, and result in a set of sample specific peaks. These
specific peaks can be used to infer sample-specific quantitation,
while the rest of the MS2 spectrum is used for identification.
Label-free MS1: extracted ion chromatograms
In label-free quantitation, the precursor peaks that match an
identified peptide are integrated of retention time and the area under
that extracted ion chromatogram is used to quantify that peptide in
that sample.
In SILAc quantitation, sample are grown in a medium that contains
heavy amino acids (typically arginine and lysine). All proteins gown
in this heavy growth medium contain the heavy form of these amino
acids. Two samples, one grown in heavy medium, and one grown in normal
(light) medium are then combined and analysed together. The heavy
peptides precursor peaks are systematically shifted compared to the
light ones, and the ratio between the height of a heavy and light
peaks can be used to calculate peptide and protein fold-changes.
Figure: credit Wikimedia Commons.
Statistical challenge These different quantitation
techniques come with their respective and distinct challenges, such as
large quantities of raw data processing, data transformation and
normalisation, and different underlying statistical models for the
quantitative data (count data for spectral counting, continuous data for
the others).
Exercise
As the name of the file implies, the data is an iTRAQ-based isobar
quantitation experiment. We can visualise the reporter peaks as
follows:
plotSpectra(spi[14])
Visualisation of the iTRAQ reporter peaks.
plotSpectra(spi[14], xlim = c(114, 117.5))
Visualisation of the iTRAQ reporter peaks.
From the closed-up figure, we see that the 4 reporter ions indicate
that there was slightly less of that peptide in the two first samples,
i.e. those labelled with the 114 and 115 reporters.
The rest of the spectrum (i.e m/z > 120) correspond to the dissociated
fragment ions, used for identification.
Quantitative data processing
Statistical challenge Quantitative data processing
come with numerous statistical challenges such as missing data handling,
aggregation of quantitation data, data normalisation, … and the effect
of these on the downstream statistical tests.
We will use the QFeatures package to handle
quantitative data.
library("QFeatures")
Let’s start by loading a PSM-level quantitation table into a
QFeatures object:
## Get columns with quantitative data
## (others contain metadata)
(quantCols <- grep("^X", colnames(hlpsms)))
## [1] 1 2 3 4 5 6 7 8 9 10
## Parse and load the data
qf <- readQFeatures(hlpsms, quantCols = quantCols,
name = "psms")
## Checking arguments.
## Loading data as a 'SummarizedExperiment' object.
## Formatting sample annotations (colData).
## Formatting data as a 'QFeatures' object.
qf
## An instance of class QFeatures (type: bulk) with 1 set:
##
## [1] psms: SummarizedExperiment with 3010 rows and 10 columns
In label-free experiments, given that each sample is acquired
independently, the proportion of missing values can be as high several
tens of percent. In such situations, removing features/rows with
missing values can be too harsh. Imputation is possible (see the
impute() function and documentation), albeit tricky, as different
mechanisms can be responsible for missing value that appear either at
random or not at random (Lazar et al. 2016).
In MS-baesd proteomics, we don’t acquire data at the protein level,
given that we cleave our peptides. We thus need to aggregate
spectrum-level quantitation values into peptide and protein-level data
using.
Below, we illustrate this aggregating values using median aggregation
and the aggregateFeatures function, first from PSMs to peptides,
then from peptides to proteins.
qf <- qf |>
aggregateFeatures(
"psms",
name = "peptides",
fcol = "Sequence",
fun = colMedians) |>
aggregateFeatures(
"peptides",
name = "proteins",
fcol = "ProteinDescriptions",
fun = colMedians)
## Your row data contain missing values. Please read the relevant
## section(s) in the aggregateFeatures manual page regarding the effects
## of missing values on data aggregation.
## Aggregated: 1/1
## Aggregated: 1/1
qf
## An instance of class QFeatures (type: bulk) with 3 sets:
##
## [1] psms: SummarizedExperiment with 3010 rows and 10 columns
## [2] peptides: SummarizedExperiment with 2923 rows and 10 columns
## [3] proteins: SummarizedExperiment with 1600 rows and 10 columns
See ?aggregateFeatures for more details.
Following on from here, many data processing such as normalisation,
non-specific filtering, and hypothesis testing is very similar to
other omics data.
You can find more details in the many Bioconductor package
vignettes. The R for Mass Spectrometry
book provide further
details and state-of-the-art infrastructure for MS and proteomics data
analysis.
Applications in statistical learning
Proteomics can be used for many applications. Simple identification
and quantitation of proteins in various different biological samples
to identify differentially abundant proteins (this is possible in bulk
and single-cells - see the scp package), but
also analysis of post-translational modifications, interactions
between proteins, and protein sub-cellular localisation, as
illustrated below.
The pRoloc package.
In this section, I illustrate a use case of mass spectrometry-based
proteomics to infer sub-cellular protein localisation and the
application of statistical learning. This section requires the
pRoloc and pRolocdata packages (L. Gatto et al. 2014; Breckels et al. 2016). Both
packages can be installed with the BiocManager::install() function,
as illustrated above.
library("pRoloc")
library("pRolocdata")
The hyperLOPIT2015 data contains localisation data for 5032 proteins
in mouse embryonic stem cells (Christoforou et al. 2016). The protein
information is collected along a 20 fraction density gradient - our
data matrix has dimensions 5032 rows and 20 columns. We use PCA to
easily visualise it in 2 dimensions.
data(hyperLOPIT2015)
## set colours
setStockcol(paste0(getStockcol(), 80))
## produce PCA plot
plot2D(hyperLOPIT2015)
## Warning in prcomp.default(x = structure(c(0.028, 0.039, 0.021, 0.12, 0.055, :
## partial argument match of 'scale' to 'scale.'
Mouse stem cell spatial proteomics data from Christoforou et al.
Each dot on the PCA plot represents a single protein. The protein is
either of unknown localisation, and represented by a grey circle, or
is of known localisation (and called an organelle marker) and is
coloured according to its expected sub-cellular localisation.
In the next figure, we have trained an support vector machine (SVM)
model and classified proteins of unknown localisation using an SVM
model. The size of the dots are scaled according to the classifier
score and only assignments corresponding to a false discovery rate of
5% have been considered.
Breckels, L M, C M Mulvey, K S Lilley, and L Gatto. 2016. “A Bioconductor Workflow for Processing and Analysing Spatial Proteomics Data.”F1000Res 5: 2926. https://doi.org/10.12688/f1000research.10411.2.
Christoforou, A, C M Mulvey, L M Breckels, A Geladaki, T Hurrell, P C Hayward, T Naake, et al. 2016. “A Draft Map of the Mouse Pluripotent Stem Cell Spatial Proteome.”Nat Commun 7 (January): 8992. https://doi.org/10.1038/ncomms9992.
Deflandre, Guillaume, Sebastian Gibb, and Laurent Gatto. 2026. “PSMatch: An R/Bioconductor Package to Explore Proteomics Identification Data.” OSF Preprints. osf.io/62v9p_v3.
Gatto, L, L M Breckels, S Wieczorek, T Burger, and K S Lilley. 2014. “Mass-Spectrometry-Based Spatial Proteomics Data Analysis Using pRoloc and pRolocdata.”Bioinformatics 30 (9): 1322–24. https://doi.org/10.1093/bioinformatics/btu013.
Gatto, L, and A Christoforou. 2014. “Using R and Bioconductor for Proteomics Data Analysis.”Biochim. Biophys. Acta 1844 (1 Pt A): 42–51.
Huber, W, V J Carey, R Gentleman, S Anders, M Carlson, B S Carvalho, H C Bravo, et al. 2015. “Orchestrating High-Throughput Genomic Analysis with Bioconductor.”Nat. Methods 12 (2): 115–21.
Käll, L, J D Storey, M J MacCoss, and W S Noble. 2008. “Posterior Error Probabilities and False Discovery Rates: Two Sides of the Same Coin.”J Proteome Res 7 (1): 40–44. https://doi.org/10.1021/pr700739d.
Lazar, C, L Gatto, M Ferro, C Bruley, and T Burger. 2016. “Accounting for the Multiple Natures of Missing Values in Label-Free Quantitative Proteomics Data Sets to Compare Imputation Strategies.”J Proteome Res 15 (4): 1116–25. https://doi.org/10.1021/acs.jproteome.5b00981.
Nesvizhskii, A I, and R Aebersold. 2005. “Interpretation of Shotgun Proteomic Data: The Protein Inference Problem.”Mol Cell Proteomics 4 (10): 1419–40. https://doi.org/10.1074/mcp.R500012-MCP200.
Rainer, Johannes, Marilyn De Graeve, Philippine Louail, Sebastian Gibb, and Laurent Gatto. 2026. “An Open Infrastructure for Mass Spectrometry Data in R,” February. osf.io/cwt2v_v3.
Sinha, Ankit, and Matthias Mann. 2020. “A beginner’s guide to mass spectrometry–based proteomics.”The Biochemist, September. https://doi.org/10.1042/BIO20200057.
Information flow in biological systems (Source Central dogma of biology on Wikipedia).
The R for Mass Spectrometry intiative.
The 
