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.

# 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

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)

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

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 …

# 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

## 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)

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 
  

Acknowledgements

This package would not have been possible without the very dedicated and hard work of my collaborator Huilin Li at Sun Yat-Sen University in China.

The author wants to acknowledge the support by the IGBMC (CNRS UMR 7104, Inserm U 1258, UdS), the proteomics platform of the IGBMC, CNRS, IGBMC, Université de Strasbourg and Inserm.

Appendix: Session-Info

#> R version 4.0.2 (2020-06-22)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 10 x64 (build 18363)
#> 
#> Matrix products: default
#> 
#> locale:
#> [1] LC_COLLATE=C                   LC_CTYPE=French_France.1252   
#> [3] LC_MONETARY=French_France.1252 LC_NUMERIC=C                  
#> [5] LC_TIME=French_France.1252    
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] wrTopDownFrag_1.0.2 wrProteo_1.1.4      wrMisc_1.3.0       
#> 
#> loaded via a namespace (and not attached):
#>  [1] compiler_4.0.2  magrittr_1.5    htmltools_0.5.0 tools_4.0.2    
#>  [5] yaml_2.2.1      stringi_1.4.6   rmarkdown_2.3   knitr_1.29     
#>  [9] stringr_1.4.0   digest_0.6.25   xfun_0.16       rlang_0.4.7    
#> [13] evaluate_0.14