Theoretic spectra generation of SIP-labeled peptide

library(Aerith)
library(ggplot2)

Overview

Stable isotope probing (SIP) experiments benefit from reliable theoretical spectra to interpret precursor and fragment ions. This vignette walks through the Aerith workflow for generating SIP-aware peptide spectra, emphasizing how each function supports isotope-resolved analysis. We start with unlabeled peptides and progress to labeled scenarios, outlining key parameters, interpretation tips, and best practices. Aerith integrates binomial-based isotope modeling, fragment intensity estimation, and visualization in a single toolkit, aligning with current SIP proteomics methodologies that rely on high accuracy precursor modeling and fragment annotation.

Compiled parameters and atom sources

Aerith uses a C++ AerithParameters object to store compiled chemistry parameters. These define real C, H, O, N, P, and S isotope distributions, residue and PTM formulas, fixed carbamidomethylation, and SIP scoring defaults. There is no configuration-file reader or file-based override. The former generateOneCFG and generateCFGs functions have been removed.

parameters <- getAerithParameters()
parameters$chemistryProfile
#> [1] "sipros5/source-aware-cam-tryptic-water/v1"
parameters$fixedPtms
#> [1] "carbamidomethyl"
parameters$residues$reagent["C", ]
#> C H O N P S 
#> 2 3 1 1 0 0

Only biosynthetic atoms follow the requested SIP abundance. Reagent atoms and atoms introduced by digestion solvent keep their natural isotope distributions. IAA’s net C2H3NO addition belongs to the reagent pool, and terminal H2O belongs to the digestion-solvent pool. Oxidation adds reagent oxygen; deamidation removes biosynthetic H/N and adds reagent oxygen. Phosphorylation uses real phosphorus. This bookkeeping does not model subsequent solvent exchange.

calPepAtomCount("C")             # Total C5 H10 O3 N2 S, including termini
#>   C  H O N P S
#> 1 5 10 3 2 0 1
calPepAtomCount("C", "sip")      # C3 H5 O N S
#>   C H O N P S
#> 1 3 5 1 1 0 1
calPepAtomCount("C", "reagent")  # IAA: C2 H3 O N
#>   C H O N P S
#> 1 2 3 1 1 0 0
calPepAtomCount("C", "solvent")  # Terminal water: H2 O
#>   C H O N P S
#> 1 0 2 1 0 0 0
calPepAtomCount("C", "natural")  # Reagent + solvent
#>   C H O N P S
#> 1 2 5 2 1 0 0

An explicit C/ is a zero-delta alias for the already blocked cysteine. The precursor and fragment calculators use the same sourced formulas; b ions contain no terminal water and y ions retain natural terminal water. Negative PTM deltas are combined with the residue before calculating its isotope envelope. Matched peak abundance estimation subtracts the expected natural background and counts only biosynthetic atoms as labelable.

Unsupported isotope names, non-finite abundances, abundances outside [0, 1], and unknown residues/PTMs raise errors. Each calculation starts from the compiled parameters, so earlier calls cannot alter the next calculation’s natural pools. Use precursor_peak_calculator_DIY to calculate a precursor spectrum from its actual sourced peptide composition.

Unlabeled Spectra

For unlabeled peptides, Aerith estimates precursor isotope distributions, calculates neutral masses, and visualizes theoretical spectra that reflect expected instrument observations.

Get precursor mass

The precursor_peak_calculator function returns the isotopic distribution for a peptide sequence, while calPepAtomCount reports elemental composition derived from the amino acid content. The calPepPrecursorMass function uses the nominal shift of the isotope envelope’s most abundant peak and its mean isotope spacing to estimate a representative precursor peak mass given a target isotope and its natural abundance. The sequence parameter accepts standard one-letter amino acid codes, and the isotope argument takes values such as "C13" or "N15". Adjust the abundance parameter to match experimental labeling levels; use natural abundance (for example 0.0107 for carbon) when no enrichment is applied.

a <- precursor_peak_calculator("PEPTIDECCCC")
head(a, 5)
#>       Mass       Prob
#> 1 1439.483 0.40252107
#> 2 1440.485 0.27774526
#> 3 1441.484 0.18623021
#> 4 1442.485 0.08384066
#> 5 1443.484 0.03375525
calPepAtomCount("PEPTIDECCCC")
#>    C  H  O  N P S
#> 1 54 85 23 15 0 4
calPepPrecursorMass("PEPTIDECCCC", "C13", 0.0107)
#> [1] 1439.483

Plot precursor’s theoretical spectra

getPrecursorSpectra constructs theoretical mass spectra for a peptide across specified charge states. The charges argument accepts individual values or integer ranges, enabling quick comparison of charge state envelopes. Interpreting the plot helps assess isotope spacing and relative intensities, which should coincide with observed MS1 signals when the peptide is unlabeled.

a <- getPrecursorSpectra("PEPTIDE", 2)
plot(a) + scale_x_continuous(breaks = seq(400, 405, by = 1)) + geom_linerange(linewidth = 0.2)
#> Scale for x is already present.
#> Adding another scale for x, which will replace the existing scale.

a <- getPrecursorSpectra("PEPTIDE", 1:2)
plot(a)

The resulting plots show the theoretical mass-to-charge distribution for the provided sequence. A narrow spacing with consistent intensity decay indicates a well-resolved isotope cluster typical for low charge state precursors.

Get peptide fragments’ mass

BYion_peak_calculator_DIY enumerates b/y fragment ions for the peptide, incorporating isotope substitution probabilities. Provide the desired isotope label and its enrichment probability to tailor the output to experimental conditions. For unlabeled peptides, set the abundance to the natural isotope frequency to obtain expected fragment masses.

df <- BYion_peak_calculator_DIY("PEPTIDE", "C13", 0.0107)
head(df, 5)
#>       Mass         Prob Kind
#> 1 147.0532 9.340350e-01   Y1
#> 2 148.0562 5.635132e-02   Y1
#> 3 149.0577 9.097310e-03   Y1
#> 4 150.0605 4.795415e-04   Y1
#> 5 151.0622 3.523521e-05   Y1

Plot peptide fragments’ theoretical spectra

getSipBYionSpectra generates labeled or unlabeled fragment spectra, while plotSipBYionLabel annotates fragment peaks. The charge state range controls which fragments appear, and the isotope count parameter governs the number of labeled atoms considered. Use higher charge states when analyzing fragmentation data collected in higher energy dissociation experiments.

a <- getSipBYionSpectra("KHRIPCDRK", "C13", 0.0107, 1:2, 2)
plot(a) + plotSipBYionLabel(a)

a <- getSipBYionSpectra("KHRIPCDRK", "C13", 0.0107, 1:2, 0)
plot(a) + plotSipBYionLabel(a)

a <- getSipBYionSpectra("KHRIPCDRK", "C13", 0.0107, 1:2, 2:3)
plot(a) + plotSipBYionLabel(a)

These plots highlight how isotope incorporation alters fragment mass distributions. Observe the position and intensity changes when varying isotope counts to anticipate how labeled fragments manifest in MS/MS spectra.

Labeled Spectra

SIP experiments enrich specific isotopes within peptides, shifting precursor and fragment masses. Aerith accommodates enrichment levels from low to complete labeling, helping quantify incorporation levels and differentiate labeled populations.

Get precursor mass

precursor_peak_calculator_DIY accepts user-defined isotope abundances, making it suitable for labeled designs. When combined with calPepPrecursorMass, users can contrast natural abundance and enriched scenarios. Selecting the isotope (Atom) and enrichment (Prob) parameters should reflect experimental labeling; for example, set Prob close to 1 for near-complete enrichment or match pulse-labeling levels (e.g., 0.55) for partial incorporation.

a <- precursor_peak_calculator_DIY("PEPTIDECCCC", "N15", 0.55)
head(a, 5)
#>       Mass         Prob
#> 1 1439.483 6.423003e-05
#> 2 1440.480 9.052470e-04
#> 3 1441.477 5.865912e-03
#> 4 1442.475 2.316484e-02
#> 5 1443.472 6.233180e-02
calPepPrecursorMass("PEPTIDECCCC", "C13", 0.55)
#> [1] 1465.567

The following calls illustrate how precursor mass estimates respond to different isotope species and enrichment levels. Interpret shifts in the resulting masses to gauge the labeling effect and to set accurate extraction windows when processing LC-MS/MS data.

calPepAtomCount("PEPTIDECCCC")
#>    C  H  O  N P S
#> 1 54 85 23 15 0 4
calPepPrecursorMass("PEPTIDECCCC", "C13", 0.0107)
#> [1] 1439.483
calPepPrecursorMass("PEPTIDECCCC", "C13", 0.5)
#> [1] 1463.561
calPepPrecursorMass("PEPTIDECCCC", "H2", 0.000115)
#> [1] 1439.483
calPepPrecursorMass("PEPTIDECCCC", "H2", 0.5)
#> [1] 1476.709
calPepPrecursorMass("PEPTIDECCCC", "N15", 0.00368)
#> [1] 1439.483
calPepPrecursorMass("PEPTIDECCCC", "N15", 0.5)
#> [1] 1445.469
calPepPrecursorMass("PEPTIDECCCC", "O18", 0.00205)
#> [1] 1439.483
calPepPrecursorMass("PEPTIDECCCC", "O18", 0.5)
#> [1] 1457.52
calPepPrecursorMass("PEPTIDECCCC", "S34", 0.0429)
#> [1] 1439.483
calPepPrecursorMass("PEPTIDECCCC", "S34", 0.5)
#> [1] 1443.477
calPepAtomCount("PEPTIDECCCC")
#>    C  H  O  N P S
#> 1 54 85 23 15 0 4
df <- precursor_peak_calculator_DIY("PEPTIDECCCC", "C13", 0.0107)
df$Mass[which.max(df$Prob)]
#> [1] 1439.483
df <- precursor_peak_calculator_DIY("PEPTIDECCCC", "C13", 0.5)
df$Mass[which.max(df$Prob)]
#> [1] 1463.561
df <- precursor_peak_calculator_DIY("PEPTIDECCCC", "O18", 0.00205)
df$Mass[which.max(df$Prob)]
#> [1] 1439.483
df <- precursor_peak_calculator_DIY("PEPTIDECCCC", "O18", 0.5)
df$Mass[which.max(df$Prob)]
#> [1] 1457.52
df <- precursor_peak_calculator_DIY("PEPTIDECCCC", "S34", 0.0429)
df$Mass[which.max(df$Prob)]
#> [1] 1439.483
df <- precursor_peak_calculator_DIY("PEPTIDECCCC", "S34", 0.5)
df$Mass[which.max(df$Prob)]
#> [1] 1443.476

Plot precursor’s theoretical spectra

getSipPrecursorSpectra visualizes isotope envelopes under enrichment. Analyze how peak intensities redistribute as the enrichment fraction increases. Wider distributions or shifted apex masses signify successful labeling, enabling downstream quantitative comparisons.

a <- getSipPrecursorSpectra("PEPTIDE", "C13", 0.55, 2)
plot(a) + scale_x_continuous(breaks = seq(400, 420, by = 1)) + geom_linerange(linewidth = 0.2)
#> Scale for x is already present.
#> Adding another scale for x, which will replace the existing scale.

a <- getSipPrecursorSpectra("PEPTIDE", "C13", 0.55, 1:2)
plot(a)

Expect the labeled spectra to shift toward higher m/z values relative to the unlabeled counterpart. Monitoring this displacement assists in verifying labeling efficiency across charge states.

Get peptide fragments’ mass

When labeling impacts fragmentation, use BYion_peak_calculator_DIY with enriched probabilities to determine the dominant fragment masses. Filtering on Kind enables investigation of particular ion series, which is useful when tracking diagnostic fragment transitions.

df <- BYion_peak_calculator_DIY("PEPTIDE", "C13", 0.55)
head(df, 5)
#>       Mass       Prob Kind
#> 1 147.0532 0.01818803   Y1
#> 2 148.0565 0.11126279   Y1
#> 3 149.0599 0.27254240   Y1
#> 4 150.0632 0.33468961   Y1
#> 5 151.0665 0.20725340   Y1
df <- BYion_peak_calculator_DIY("PEPTIDECCCC", "C13", 0.0107)
df <- df[which(df$Kind == "Y6"), ]
df$Mass[which.max(df$Prob)]
#> [1] 902.2027
df <- BYion_peak_calculator_DIY("PEPTIDECCCC", "C13", 0.5)
df <- df[which(df$Kind == "Y6"), ]
df$Mass[which.max(df$Prob)]
#> [1] 913.2378
df <- BYion_peak_calculator_DIY("PEPTIDECCCC", "H2", 0.000115)
df <- df[which(df$Kind == "Y6"), ]
df$Mass[which.max(df$Prob)]
#> [1] 902.2027
df <- BYion_peak_calculator_DIY("PEPTIDECCCC", "H2", 0.5)
df <- df[which(df$Kind == "Y6"), ]
df$Mass[which.max(df$Prob)]
#> [1] 919.3052
df <- BYion_peak_calculator_DIY("PEPTIDECCCC", "N15", 0.00368)
df <- df[which(df$Kind == "Y6"), ]
df$Mass[which.max(df$Prob)]
#> [1] 902.2027
df <- BYion_peak_calculator_DIY("PEPTIDECCCC", "N15", 0.5)
df <- df[which(df$Kind == "Y6"), ]
df$Mass[which.max(df$Prob)]
#> [1] 905.1953

Plot peptide fragments’ theoretical spectra

The fragment spectra for labeled peptides reveal how isotope incorporation propagates across fragment ions. Examine the annotations produced by plotSipBYionLabel to confirm which ion series carry the label and to identify shifts that support peptide identification and quantification.

a <- getSipBYionSpectra("KHRIPCDRK", "C13", 0.55, 1:2, 2)
plot(a) + plotSipBYionLabel(a)

a <- getSipBYionSpectra("KHRIPCDRK", "C13", 0.55, 1:2, 0)
plot(a) + plotSipBYionLabel(a)

a <- getSipBYionSpectra("KHRIPCDRK", "C13", 0.55, 1:2, 2:3)
plot(a) + plotSipBYionLabel(a)

Through these visualizations, Aerith showcases its advantage in modeling SIP-specific fragmentation patterns in line with current state-of-the-art approaches that couple theoretical spectra with high-resolution MS/MS data analysis.

sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] tidyr_1.3.2    ggplot2_4.0.3  stringr_1.6.0  dplyr_1.2.1    Aerith_1.1.3  
#> [6] rmarkdown_2.32
#> 
#> loaded via a namespace (and not attached):
#>  [1] sass_0.4.10          generics_0.1.4       stringi_1.8.9       
#>  [4] digest_0.6.39        magrittr_2.0.5       evaluate_1.0.5      
#>  [7] grid_4.6.1           RColorBrewer_1.1-3   fastmap_1.2.0       
#> [10] jsonlite_2.0.0       ggrepel_0.9.8        ProtGenerics_1.45.0 
#> [13] mzR_2.47.1           purrr_1.2.2          scales_1.4.0        
#> [16] codetools_0.2-20     jquerylib_0.1.4      cli_3.6.6           
#> [19] rlang_1.3.0          Biobase_2.73.2       withr_3.0.3         
#> [22] cachem_1.1.0         yaml_2.3.12          otel_0.2.0          
#> [25] tools_4.6.1          ncdf4_1.24           BiocGenerics_0.59.12
#> [28] buildtools_1.0.0     vctrs_0.7.3          R6_2.6.1            
#> [31] lifecycle_1.0.5      pkgconfig_2.0.3      pillar_1.11.1       
#> [34] bslib_0.12.0         gtable_0.3.6         glue_1.8.1          
#> [37] data.table_1.18.6.1  Rcpp_1.1.2           xfun_0.61           
#> [40] tibble_3.3.1         tidyselect_1.2.1     sys_3.4.3           
#> [43] knitr_1.52           farver_2.1.2         htmltools_0.5.9     
#> [46] maketools_1.3.2      labeling_0.4.3       compiler_4.6.1      
#> [49] S7_0.2.2