Getting started with wrTopDownFrag

Introduction

This package contains tools for the use in TopDown Proteomics. Proteomics is referred to the technique of idenifying all proteins in a given sample. Typically this is done by using mass spectrometry. This technqiue returns primarily ‘m/z’ (mass over charge) measures, in most cases this can be resolved to molecular masses.

Since mass spectrometry is rather well suited to identifying small molecules, it has become common practice to first digest proteins into smaller units (ie peptides), which can be easily identified by mass spetrometry. This technique is referred to as “bottom-up” proteomics. Further high energy fragmentation of proteins or peptides directly within the mass spetrometer also helps very much improving the identification rate (MS-MS or MS2).

To overcome some of the drawbacks associated with the bottom-up approach, more recent developments of mass spectrometers allow the identification of full length proteins (from samples with few proteins). This approach is called “top-down proteomics” (see also Chen et al, 2018, Skinner et al, 2018 or Li et al, 2018). In this context high energy random fragmentation of proteins/peptides within the mass spectrometer plays an important role, too. This approach produces fragments conainting on of the original start/end-sites (terminal fragments) and, depending on the energy settings, furher internal fragments. Of course, larger parent proteins/peptides will give even more complex patterns of internal fragments. The pattern of resulting fragments for a given precursor protein/peptide allows better identification and provides further valuable information about the (3-dimensional) conformation of the initial proteins (see also Haverland et al, 2017).

This project got started to help analyzing internal fragments from FT-ICR mass-spectrometry. At the time of beginning none or only very limited tools were available for this task (this situation is changing since 2019). Since there is already software available for transforming initial lists of m/z values into monoisotopic values, the aim was to continue further to the identification of m/z peaks after a deconvolution step to assign the most likely peptide/proptein sequence. Please refer eg to Wikipedia: monoisotopic mass for details on molecular mass and deconvolution. Initial developments were performed based on data from FT-ICR mass-spectrometry, but the overal concept may be applied to any kind of mass-spectrometry data.

When developing this package particular attention was brought to the fact that entire proteins may very well still carry embedded ions in catalytic cores or other rare modifications. With the tools pesented here, the link between parent- and children fragments was not further taken into account since we did not expect 100 percent pure parent species entering the fragmentation step.

In summary, this package aims to provide tools for the identification of proteins from monoisotopic m/z lists, and in particular, to consider and identify all possible internal fragments resulting from fragmentation performed during mass-spectromery analysis.

To get started, we need to load the packages wrMisc and wrProteo, available from CRAN. And of course we need to charge this package.

library(wrMisc)
library(wrProteo)
library(wrTopDownFrag)

Further information about the package versions used for making this vignette can be found in the appendix ‘Session-info’.

For manipulating peptide/protein sequences we will use functions for working with one-letter code amino-acid sequences provided by package wrProteo.

Nomenclature

This describes how chemical modifications on amino-acids (like oxygenation) are abbreviated and what exact chemical modification it referrs to. In term of nomenclature we’ll stick to these abbreviations :

## common settings
str(AAfragSettings())
#> List of 4
#>  $ knownMods        :List of 9
#>   ..$ noMod   : chr ""
#>   ..$ Nterm   : chr [1:3] "a" "b" "c"
#>   ..$ Cterm   : chr [1:3] "x" "y" "z"
#>   ..$ NCterm  : chr [1:2] "d" "i"
#>   ..$ any     : NULL
#>   ..$ intern  : chr [1:8] "a" "b" "c" "f" ...
#>   ..$ spcAA   : chr [1:9] "p" "h" "k" "o" ...
#>   ..$ spcNterm: chr [1:2] "g" "n"
#>   ..$ spcCterm: chr "u"
#>  $ specAAMod        :List of 9
#>   ..$ p: chr [1:3] "T" "S" "Y"
#>   ..$ h: chr [1:4] "S" "T" "E" "D"
#>   ..$ k: chr [1:4] "Q" "K" "R" "N"
#>   ..$ o: chr "M"
#>   ..$ r: chr "C"
#>   ..$ s: chr "S"
#>   ..$ m: chr [1:2] "R" "K"
#>   ..$ n: chr "Q"
#>   ..$ u: chr "R"
#>  $ modChem          : chr [1:22, 1:2] "" "a" "b" "c" ...
#>   ..- attr(*, "dimnames")=List of 2
#>   .. ..$ : NULL
#>   .. ..$ : chr [1:2] "ty" "mod"
#>  $ neutralLossOrGain: chr [1:7] "b" "d" "f" "h" ...

Here we can see that ‘a’,‘b’ and ‘c’-ions are grouped as ‘Nterm’ or that the phosporylation modification ‘p’ is taken into consideration at ‘S’,‘T’ or ‘Y’ amino-acid residues. The section $modChem describes/defines exactely how many molecules will be added or removed with the various modifications available in this package.

Molecular (mono-isotopic) masses of the atomes used here are taken the package wrProteo, intially they were taken from Unimod. They can be easily updated, if in the future, (mono-isotopic) molecular masses will be determined with higher precision (ie more digits).

Obtaining Molecular Mass For Chemical Structures (using wrProteo)

The package wrProteo is used to convert summed chemical formulas into molecular mass.

# Standard way to obtain the (monoisotopic) mass of water (H20) or a phosphorylation (PO3)
# Note that the number a molecule appears must be written in front of the molecule (no number means one occurance)
massDeFormula(c("2HO", "P3O"))
#>     2H1O     1P3O 
#> 18.01056 78.95851

# Undereith this runs (for H20):
2*.atomicMasses()["H","mono"] +.atomicMasses()["O","mono"]   # H2O
#> [1] 18.01056

Atomic masses can be calulated either as ‘average mass’ or ‘monoisotopic mass’ (Wikipedia: monoisotopic mass), the latter is commonly used in mass-spectrometry and will be used by default in this package.

Molecular mass of peptides and proteins

At this level we can compute the expected mass of uncharged proteins/peptides as defined by their size. Of couse, protein isomers (ie same total composition but different sequence) get the same mass.

# Let's define two small amino-acid sequences 
protP <- c(pepK="KPEPTIDRPEP", pepC="CEPEPTRT", pepC2="PECEPTRT")

# The sequence converted to mass  
convAASeq2mass(protP)   
#> KPEPTIDRPEP    CEPEPTRT    PECEPTRT 
#>   1277.6616    931.4069    931.4069

As mentinned, this package assumes that experimental values have already been deconvoluted, ie that only the mono-charged peak provided and isotopic patterns have been reduced to the main representative isotope. In line with this assumption, default predictions are mono-isotopic masses.

Fragmenting a Peptide/Protein -Sequence

With ‘Fragmentation’ techniques we refer to technqiues like shooting electrons or IR-waves allowing to break larger molecules into smaller molecules (of different composition). In order to check if fragmentation yields random cleavage or raher directed distribution of cleavage sites, it is necessary to predict all possible cleavage sites.

The complexity of this simple task increases about exponentially with protein size. At this level, the complexity increases so much, that only a few (longer) full length proteins can be treated at once within a reasonable amount of time. With large proteins (more than 600 aa length) this may consume considerable amounts of RAM. When designing this package care has been taken to run infractructure intensive as efficent as possible, but working with complex samples/proteomes it is still beyond technical limits.

Here a very simplified view on the theoretical number of terminal and internal fragments. Fixed modifications do not change the number of expected fragments Note, that this simplification does not include variable modifications (see also the next section).

marks <- data.frame(name=c("Ubiquitin\n76aa","Glutamate dehydrogenase 1\n501aa"),length=c(76,501))
layout(matrix(1:2,ncol=2))
plotNTheor(x=20:750, log="", mark=marks)
plotNTheor(x=20:750, log="xy", mark=marks)
mtext("log/log scale", cex=0.8, line=0.1)

Fragmenting a protein sequence

Random cleavege for a sample collection of proteins can be obtained using the functions or

protP <- c(pepK="KPEPTIDRPEP", pepC="CEPEPTRT")
## Basic output 
fragmentSeq(protP[1], minSize=3, internFragments=TRUE, pref="pepK")
#>  -> fragmentSeq : 1 out of 45 fragments not unique
#>  full.pepK.1-11   Nter.pepK.1-3   Nter.pepK.1-4   Nter.pepK.1-5   Nter.pepK.1-6 
#>   "KPEPTIDRPEP"           "KPE"          "KPEP"         "KPEPT"        "KPEPTI" 
#>   Nter.pepK.1-7   Nter.pepK.1-8   Nter.pepK.1-9  Nter.pepK.1-10  Cter.pepK.2-11 
#>       "KPEPTID"      "KPEPTIDR"     "KPEPTIDRP"    "KPEPTIDRPE"    "PEPTIDRPEP" 
#>  Cter.pepK.3-11  Cter.pepK.4-11  Cter.pepK.5-11  Cter.pepK.6-11  Cter.pepK.7-11 
#>     "EPTIDRPEP"      "PTIDRPEP"       "TIDRPEP"        "IDRPEP"         "DRPEP" 
#>  Cter.pepK.8-11  Cter.pepK.9-11 inter.pepK.2-10  inter.pepK.2-4  inter.pepK.2-5 
#>          "RPEP"           "PEP"     "PEPTIDRPE"           "PEP"          "PEPT" 
#>  inter.pepK.2-6  inter.pepK.2-7  inter.pepK.2-8  inter.pepK.2-9 inter.pepK.3-10 
#>         "PEPTI"        "PEPTID"       "PEPTIDR"      "PEPTIDRP"      "EPTIDRPE" 
#> inter.pepK.4-10 inter.pepK.5-10 inter.pepK.6-10 inter.pepK.7-10 inter.pepK.8-10 
#>       "PTIDRPE"        "TIDRPE"         "IDRPE"          "DRPE"           "RPE" 
#>  inter.pepK.3-9  inter.pepK.3-5  inter.pepK.3-6  inter.pepK.3-7  inter.pepK.3-8 
#>       "EPTIDRP"           "EPT"          "EPTI"         "EPTID"        "EPTIDR" 
#>  inter.pepK.4-9  inter.pepK.5-9  inter.pepK.6-9  inter.pepK.7-9  inter.pepK.4-8 
#>        "PTIDRP"         "TIDRP"          "IDRP"           "DRP"         "PTIDR" 
#>  inter.pepK.4-6  inter.pepK.4-7  inter.pepK.5-8  inter.pepK.6-8  inter.pepK.5-7 
#>           "PTI"          "PTID"          "TIDR"           "IDR"           "TID"

## Elaborate output
protP2 <- cbind(na=names(protP), se=protP, ma=wrProteo::convAASeq2mass(protP,seqName=TRUE))
pepT1 <- makeFragments(protTab=protP2, minFragSize=3, maxFragSize=9, internFra=TRUE)
#> -> makeFragments -> fragmentSeq : 1 out of 42 fragments not unique
#>  -> makeFragments :  due to preferential fragmentation sites discard 2 fragments,  61 remain

head(pepT1)
#>      no   seq   orig          origNa ty      seqNa       beg end  precAA tailAA
#> [1,] "49" "PTI" "KPEPTIDRPEP" "pepK" "inter" "pepK.4-6"  "4" "6"  "E"    "D"   
#> [2,] "20" "PEP" "KPEPTIDRPEP" "pepK" "Cter"  "pepK.9-11" "9" "11" "R"    NA    
#> [3,] "27" "PEP" "KPEPTIDRPEP" "pepK" "inter" "pepK.2-4"  "2" "4"  "K"    "T"   
#> [4,] "62" "PEP" "CEPEPTRT"    "pepC" "inter" "pepC.3-5"  "3" "5"  "E"    "T"   
#> [5,] "40" "EPT" "KPEPTIDRPEP" "pepK" "inter" "pepK.3-5"  "3" "5"  "P"    "I"   
#> [6,] "63" "EPT" "CEPEPTRT"    "pepC" "inter" "pepC.4-6"  "4" "6"  "P"    "R"   
#>      ambig          mass               modSpec
#> [1,] NA             "329.194522412891" ""     
#> [2,] "duplSequence" "341.158136906591" ""     
#> [3,] "duplSequence" "341.158136906591" ""     
#> [4,] "duplSequence" "341.158136906591" ""     
#> [5,] "duplSequence" "345.153051528691" ""     
#> [6,] "duplSequence" "345.153051528691" ""

dim(pepT1)
#> [1] 61 13

## The repartition between types of fragments :
table(pepT1[,"ty"])
#> 
#>  Cter  Nter  full inter 
#>    12    12     1    36

## Types of ambiguities encountered
table(pepT1[,"ambig"])
#> 
#> duplSequence      isoMass 
#>            7           10

Even such a small example gives already 61 possible peptides (without counting modfications). Of course, it is quite common to obatin a high degree of ambiguities with short peptides since they are less likely unqique.

Fixed And Variable Modifications

In real-world biology protein modifcations are common. Conceptually one can distuinguish two cases : With ‘fixed modificatons’ it is presumed that all protein moleculs of a given species (ie sequence) do carry exactely the same modification(s). Alteratively one may suggest that only a portion of the molecules for a given protein-species carry a given modification.

Fixed modifications do not increase the search space, since a given value corresponding to the change of mass gets added or subtracted. In the case of varaible modifications this increases the search space since the modificed as well as the unmodified mass will be considered when comparig to experimental masses. In particular, when multiple amino-acids on the same protein may get modified alternatively this increases the search space. Then one may consider multiple modifications, for example 0 to 2 out of 5 Serine residues in the protein sequence may carry a phosporylation …

In order to reduce the apparent complexity, several conceptual compromises have been taken. ) The presence of charged amino-acids has been used to dismiss all fragments not containing a charged amino-acid, default has been set to positive charge (K,H and R). ) When identifying variable modifications, the non-modified isoform must be identified first to take the variable modification in account.

Basic Identification Including Fixed Modifications

Now we are ready to compare a list of experimental m/z values to a set of protein sequences …

Tolerance

In mass spectrometry it is common to use relative tolerance-limits when it comes to identification as ‘ppm’. This means that peaks not having any predicted peptides/ions in a predefined ppm-range will be omitted. When the search space gets very crowded, ie when many peptide-fragments are predicted, there is of course a considerable risk that some predicted fragments/ions are so close that multiple predicted peptides/ions have to be consiered in a a given ppm-range.

Identification Example

First let’s make a little toy example using a hypothetic protein/peptide sequence. The molecular masses below have been predicted using the Prospector tool from UCSF. We’ll add than a litle bit of random noise as a simultation for experimental data. The table below is not exhaustive for all potentially occurring fragments, as experimental data would not be expected to be complete either. As modifications to test some cases of loss of water (H20) and loss of ammonia (NH3) were prepared.

# The toy peptide/protein sequnce
protP <- c(pepK="KPEPTIDRPEP", pepC="CEPEPTRT")

obsMass1 <- cbind(a=c(424.2554,525.3031,638.3872,753.4141,909.5152,1006.5680,1135.6106),
  b=c(452.2504,553.2980,666.3821,781.4090,937.5102,1034.5629,1163.6055),
  x=c(524.2463,639.2733,752.3573,853.4050,950.4578,1079.5004,1176.5531),
  y=c(498.2671,613.2940,726.3781,827.4258,924.4785,1053.5211,1150.5739),
  bdH=c(434.2398,535.2875,648.3715,763.3985,919.4996,1016.5524,1145.5949),
  ydH=c(480.2565,1132.5633,595.2835,708.3675,809.4152,906.4680,1035.5106),
  bi=c(498.2307,583.3198,583.3198,611.3148,680.3726,712.3624,712.3624),
  bidH=c(662.3620,694.3519,694.3519,791.4046,791.4046,791.4046,888.4574),
  bidN=c(663.3461,695.3359,695.3359,792.3886,792.3886,792.3886,889.4414),
  ai=c(652.3777,684.3675,684.3675,781.4203,781.4203,781.4203,878.4730) )
rownames(obsMass1) <- c("P","T","I","D","R","P","E")      # only for N-term (a & b)

## This example contains several iso-mass cases
length(obsMass1)
#> [1] 70
length(unique(as.numeric(obsMass1)))
#> [1] 59

Iso-peptides will show up with the same mass :

## have same mass   
##       480.2201   DRPE-H2O & y4-H2O
##       583.3198   PTIDR & TIDRP ;         565.3093    PTIDR-H2O & TIDRP-H2O 
##       809.4152   y7-H2O  & EPTIDRP  
##       684.3675   TIDRPE-CO & EPTIDR-CO ; 694.3519    EPTIDR-H2O & TIDRPE-H2O
##       712.3624   TIDRPE & EPTIDR

# Now we'll add some random noise
set.seed(2020)
obsMass2 <- as.numeric(obsMass1)
obsMass2 <- obsMass2 + runif(length(obsMass2), min=-2e-4, max=2e-4)

Now let’s use these numbers as experimental values and compare them to all theoretical values : For example, the peptide-sequences ‘PTIDR’ and ‘TIDRP’ may be derived from our first toy-sequence and contain exactely the same atoms and thus will have iso-masses. Without any further information it is impossible to know which one of them is/was truly present in a given sample. Thus, the corresponding mass will be called an ambiguous identification.


identP1 <- identifFixedModif(prot=protP[1], expMass=obsMass2, minFragSize=3, 
  maxFragSize=7, modTy=list(basMod=c("b","y")))     # should find 10term +10inter
#> -> identifFixedModif -> makeFragments -> fragmentSeq : 1 out of 35 fragments not unique
#> -> identifFixedModif -> makeFragments :  due to preferential fragmentation sites discard 2 fragments,  33 remain
#>  -> identifFixedModif :  removing 11 out of 33 (initial) peptides not containing any charge-catching residues
#>  -> identifFixedModif : 48 out of 70 experim masses in range of 22 predicted (max 828)
#> -> identifFixedModif -> findCloseMatch :  note that some names do overlap (beware not to confuse...) !
#>  -> identifFixedModif :  1st pass: compare 22 predicted (incl 4 ambiguous : isoMass
#>      with  48 input (measured) masses, found  16 groups of matches to experimental masses (total=19))
 
## This function returns a list   
str(identP1)  
#> List of 6
#>  $ massMatch  :List of 16
#>   ..$ 2 : Named num 0.00668
#>   .. ..- attr(*, "names")= chr "8"
#>   ..$ 3 : Named num 0.43
#>   .. ..- attr(*, "names")= chr "9"
#>   ..$ 4 : Named num -0.0681
#>   .. ..- attr(*, "names")= chr "10"
#>   ..$ 5 : Named num -0.0768
#>   .. ..- attr(*, "names")= chr "11"
#>   ..$ 9 : Named num 0.285
#>   .. ..- attr(*, "names")= chr "22"
#>   ..$ 19: Named num -0.317
#>   .. ..- attr(*, "names")= chr "43"
#>   ..$ 30: Named num [1:2] -0.2007 -0.0225
#>   .. ..- attr(*, "names")= chr [1:2] "44" "45"
#>   ..$ 8 : Named num -0.182
#>   .. ..- attr(*, "names")= chr "23"
#>   ..$ 18: Named num -0.126
#>   .. ..- attr(*, "names")= chr "46"
#>   ..$ 26: Named num 0.0181
#>   .. ..- attr(*, "names")= chr "47"
#>   ..$ 7 : Named num -0.271
#>   .. ..- attr(*, "names")= chr "24"
#>   ..$ 17: Named num [1:2] -0.0707 -0.1253
#>   .. ..- attr(*, "names")= chr [1:2] "48" "49"
#>   ..$ 25: Named num [1:2] -0.0707 -0.1253
#>   .. ..- attr(*, "names")= chr [1:2] "48" "49"
#>   ..$ 6 : Named num -0.236
#>   .. ..- attr(*, "names")= chr "25"
#>   ..$ 15: Named num 0.13
#>   .. ..- attr(*, "names")= chr "40"
#>   ..$ 16: Named num 0.13
#>   .. ..- attr(*, "names")= chr "40"
#>  $ preMa      : chr [1:22, 1:14] "1" "2" "3" "4" ...
#>   ..- attr(*, "dimnames")=List of 2
#>   .. ..$ : chr [1:22] "pepK.1-3.b" "pepK.1-4.b" "pepK.1-5.b" "pepK.1-6.b" ...
#>   .. ..$ : chr [1:14] "no" "seq" "orig" "origNa" ...
#>  $ pepTab     : chr [1:22, 1:13] "1" "2" "3" "4" ...
#>   ..- attr(*, "dimnames")=List of 2
#>   .. ..$ : NULL
#>   .. ..$ : chr [1:13] "no" "seq" "orig" "origNa" ...
#>  $ recalibFact: num 0
#>  $ recalibData: NULL
#>  $ docTi      : Named num [1:7] 1.6e+09 1.6e+09 1.6e+09 1.6e+09 1.6e+09 ...
#>   ..- attr(*, "names")= chr [1:7] "ini_identifFixedModif" "makeFragments" "countPotModifAAs" "addMassModif" ...
## $masMatch identifies each

However, due to real-world imprecision during the process of measuring m/z, here we used a 5ppm default tolerance.

This brings us to one of the reasons why fragmentation is so important in proteomics : Without fragmentation mass spectrometry of entire proteins is in big trouble to resolve most of naturally occuring iso-variants or even related proteins.

Due to fragmentation numerous/may overlapping fragments will occur and will finally allow reconstructing the most parts of original protein.

Thank you for you interest in this package. This package is still under development, new functions will be added to the next version.