Using the Python API

This tutorial shows how to use the MS²Rescore Python API for each step of the rescoring process individually. This is useful if you want to customize rescoring for your own Python workflow or if you want to understand how MS²Rescore works. Note that the full MS²Rescore workflow is also available from Python with the single function call ms2rescore.rescore().

[1]:

import logging
import plotly.io

logging.basicConfig(level=logging.INFO)
plotly.io.renderers.default = "plotly_mimetype+notebook"

Reading and parsing peptide-spectrum matches

[2]:

from psm_utils.io import read_file

from ms2rescore.report.charts import score_histogram

Reading the PSM file

MS²Rescore is fully centered around the use of a psm_utils PSMList. This is a unified data representation of PSMs and their various attributes. Internally, it is simply a list of Pydantic data classes which represent PSMs. With the submodule psm_utils.io, we can read PSMs from a variety of file formats. Here, we will read a PSM file in the MaxQuant msms.txt format.

Importantly, for rescoring, the PSM file must contain all target and decoy PSMs, including PSMs that did not pass the FDR threshold. Most search engines must be specifically configured to return all PSMs without FDR filtering.

[3]:

psm_list = read_file("msms.txt", filetype="msms")
psm_list.calculate_qvalues()

For a quick inspection, we can format the PSM list as a Pandas dataframe and display the first few rows:

[4]:

psm_list.to_dataframe().head()

[4]:
peptidoform spectrum_id run collection spectrum is_decoy score qvalue pep precursor_mz retention_time ion_mobility protein_list rank source provenance_data metadata rescoring_features
0 ((A, None), (A, None), (A, None), (A, None), (... 4703 20161213_NGHF_DBJ_SA_Exp3A_HeLa_1ug_7min_15000_02 None None False 107.660 0.001566 0.001517 478.77982 5.2007 None [P36578, H3BM89, H3BU31] None msms {'msms_filename': 'msms.txt'} {'Scan index': '3698', 'Sequence': 'AAAAAAALQA... {}
1 ((A, None), (A, None), (A, None), (A, None), (... 13572 20161213_NGHF_DBJ_SA_Exp3A_HeLa_1ug_7min_15000_02 None None False 107.740 0.001566 0.004931 915.15197 11.8470 None [O00410, E7ETV3, E7EQT5, C9JZD8] None msms {'msms_filename': 'msms.txt'} {'Scan index': '11885', 'Sequence': 'AAAAAEQQQ... {}
2 ((A, None), (A, None), (A, None), (A, None), (... 13366 20161213_NGHF_DBJ_SA_Exp3A_HeLa_1ug_7min_15000_02 None None False 137.890 0.000673 0.000493 957.85029 11.6900 None [O00410, E7ETV3, E7EQT5, C9JZD8] None msms {'msms_filename': 'msms.txt'} {'Scan index': '11695', 'Sequence': 'AAAAAEQQQ... {}
3 ((A, None), (A, None), (A, None), (A, None), (... 505 20161213_NGHF_DBJ_SA_Exp3A_HeLa_1ug_7min_15000_02 None None False 22.641 0.487532 0.142020 585.28653 0.5178 None [E9PJF0, E9PQW4, P27361] None msms {'msms_filename': 'msms.txt'} {'Scan index': '419', 'Sequence': 'AAAAAQGGGGG... {}
4 ((A, None), (A, None), (A, None), (A, None), (... 6589 20161213_NGHF_DBJ_SA_Exp3A_HeLa_1ug_7min_15000_02 None None False 89.403 0.003778 0.046504 823.87389 6.6105 None [Q9P258] None msms {'msms_filename': 'msms.txt'} {'Scan index': '5439', 'Sequence': 'AAAAAWEEPS... {}

We can also directly plot the current PSM score distributions:

[5]:

score_histogram(psm_list)

Parsing modification names

While psm_utils could take care of all file parsing, we must still map the amino acid modification names that were used by the search engine to ones that are recognized by tools such as MS²PIP and DeepLC. This includes:

  • Names as used in the Unimod or PSI-MOD databases

  • Accession numbers as used in the Unimod or PSI-MOD databases

  • Chemical formulas

  • Mass shifts in Da

Note that, for instance DeepLC, requires a chemical formula to encode modifications. It will not be able to correctly encode modifications if only the mass shift is provided. It is therefore always preferred to provide a name/accession of a database where this information can be retrieved, or provide the chemical formula directly.

If a chemical formula is provided, other tools, such as MS²PIP, can use it to derive the correct mass shift.

To map modification names, simply provide a dictionary to the psm_list.rename_modifications method.

[6]:

psm_list.rename_modifications({
    "gl": "U:Gln->pyro-Glu",
    "ox": "U:Oxidation",
    "ac": "U:Acetylation",
    "de": "U:Deamidation",
})

Assigning fixed modifications

Some search engines (although not many) do not report fixed (sometimes also called static) modifications in their output PSM files. This is for instance the case for MaxQuant. With the method psm_list.add_fixed_modifications, we can systematically assign fixed modifications to their amino acid targets, as was done during the search.

Note that add_fixed_modifications adds the modification in the ProForma 2.0 encoding as a prefix. To fully apply the modifications in the sequence, use the psm_list.apply_fixed_modifications method. A PSM ACDE/2 would therefore go from <[U:4]@C>ACDE/2 to AC[U:4]DE/2.

[7]:

psm_list.add_fixed_modifications([("U:Carbamidomethyl", ["C"])])
psm_list.apply_fixed_modifications()

Adding rescoring features

[8]:

import pandas as pd

from ms2rescore.feature_generators.basic import BasicFeatureGenerator
from ms2rescore.feature_generators.ms2 import MS2FeatureGenerator
from ms2rescore.feature_generators.ms2pip import MS2PIPFeatureGenerator
from ms2rescore.feature_generators.deeplc import DeepLCFeatureGenerator

The PSM list, as read from a MaxQuant msms.txt file, currently does not contain any rescoring features:

[9]:

pd.DataFrame(list(psm_list["rescoring_features"]))

[9]:
0
1
2
3
4
...
35551
35552
35553
35554
35555

35556 rows × 0 columns

Note that psm_list[rescoring_features] returns a Numpy array of dictionaries, but we can use Pandas to format it as a table.

Getting precursor information and annotating MS2 spectra

Several feature generators require data from the spectrum files, such as retention times, precursor m/z values, or the MS2 spectra themselves. Rather than having each feature generator parse and annotate spectra independently, MS²Rescore v4.0+ loads and annotates the required data once upfront:

  1. add_precursor_values reads spectrum files and attaches raw MS2Spectrum objects and precursor metadata (RT, precursor m/z) to each PSM.

  2. annotate_spectra then annotates each spectrum with its PSM’s peptidoform, replacing the raw spectrum with an AnnotatedMS2Spectrum. All subsequent feature generators (MS²PIP, MS2, etc.) can reuse these annotations directly, so annotation only happens once.

The following two parameters are especially important:

  • spectrum_path: Path to the original spectrum files. This can be a single file, or a directory containing multiple files. Supported formats: MGF, mzML.

  • spectrum_id_pattern: A regular expression pattern to match the spectrum ID from the PSM file to the spectrum identifier in the spectrum file.

A key aspect of spectrum parsing is matching peptide identifications to the correct spectra. Many search engines alter spectrum identifiers in their output files, so the identifiers in the PSM file may not directly match those in the spectrum file. For example, MaxQuant reports scan numbers rather than full spectrum titles. If the spectrum files were converted to MGF with ThermoRawFileParser, we can extract these scan numbers from the spectrum TITLE entries using a regular expression pattern. The following pattern extracts the scan number from this example spectrum title:

spectrum_id_pattern = r"scan=(\d+)"

mzspec=20161213_NGHF_DBJ_SA_Exp3A_HeLa_1ug_7min_15000_02.raw: controllerType=0 controllerNumber=1 scan=2

We can test this with the first spectrum in the file:

[10]:

from pyteomics.mgf import read as read_mgf

for spectrum in read_mgf("20161213_NGHF_DBJ_SA_Exp3A_HeLa_1ug_7min_15000_02.mgf"):
    print(spectrum["params"]["title"])  # noqa T201
    break

mzspec=20161213_NGHF_DBJ_SA_Exp3A_HeLa_1ug_7min_15000_02.raw: controllerType=0 controllerNumber=1 scan=2

[11]:

import re
spectrum_id = re.match(r".*scan=(\d+)$", spectrum["params"]["title"]).group(1)
print(spectrum_id)  # noqa T201

2

We can use this scan number to retrieve its corresponding PSM:

[12]:

example_psm = psm_list[psm_list["spectrum_id"] == spectrum_id][0]
example_psm

[12]:

PSM(peptidoform=Peptidoform('ELMKLAPR/2'), spectrum_id='2', run='20161213_NGHF_DBJ_SA_Exp3A_HeLa_1ug_7min_15000_02', collection=None, spectrum=None, is_decoy=False, score=58.194, qvalue=np.float64(0.04256930768065444), pep=0.11486, precursor_mz=479.28114, retention_time=0.0054151, ion_mobility=None, protein_list=['K7ER25', 'K7EJY1', 'Q7Z3S7', 'E7EUE0'], rank=None, source='msms', provenance_data={'msms_filename': 'msms.txt'}, metadata={'Scan index': '0', 'Sequence': 'ELMKLAPR', 'Length': '8', 'Missed cleavages': '1', 'Modifications': 'Unmodified', 'Deamidation (NQ) Probabilities': '', 'Oxidation (M) Probabilities': '', 'Deamidation (NQ) Score Diffs': '', 'Oxidation (M) Score Diffs': '', 'Acetyl (Protein N-term)': '0', 'Deamidation (NQ)': '0', 'Gln->pyro-Glu': '0', 'Oxidation (M)': '0', 'Gene Names': 'CACNA2D4', 'Protein Names': 'Voltage-dependent calcium channel subunit alpha-2/delta-4;Voltage-dependent calcium channel subunit alpha-2-4;Voltage-dependent calcium channel subunit delta-4', 'Fragmentation': 'HCD', 'Mass analyzer': 'FTMS', 'Type': 'MSMS', 'Scan event number': '1', 'Isotope index': '', 'Mass': '956.54772', 'Mass error [ppm]': 'NaN', 'Mass error [Da]': 'NaN', 'Simple mass error [ppm]': '10020052000', 'Delta score': '12.312', 'Score diff': 'NaN', 'Localization prob': 'NaN', 'Combinatorics': '1', 'PIF': '0', 'Fraction of total spectrum': '0', 'Base peak fraction': '0', 'Precursor Full ScanNumber': '-1', 'Precursor Intensity': 'NaN', 'Precursor Apex Fraction': 'NaN', 'Precursor Apex Offset': '0', 'Precursor Apex Offset Time': 'NaN', 'Matches': 'y1;y2;y3;y1-NH3;y2-NH3;y3-NH3;b5;b6;b7;b5-H2O;b6-H2O', 'Intensities': '1301231.4;1312401.4;842082.2;527518.5;3855320.6;670603.1;964535.9;2367717.9;197335.6;30821.9;199087.7', 'Mass Deviations [Da]': '-0.000415732784006195;-0.000550409580171163;-0.00104542525048146;-0.000307003389082183;-0.000491550706556154;-0.000489022332658351;0.00622885583993593;0.00615890510266581;0.00973991597902568;0.0034132017480033;0.00802109195342382', 'Mass Deviations [ppm]': '-2.37399659991276;-2.02228385487352;-3.04602322334689;-1.94191995899886;-1.92654936588088;-1.49922797285001;10.1225059202731;8.97296777490114;12.4323442241631;5.71400675006929;12.0009400088042', 'Masses': '175.119367913784;272.17226644258;343.20987524605;158.092710082889;255.145658482207;326.182769741633;615.34721629166;686.384400030197;783.433582871321;597.339467259452;668.371973157047', 'Number of Matches': '11', 'Intensity coverage': '0.0871847693164245', 'Peak coverage': '0.0728476821192053', 'Neutral loss level': 'None', 'ETD identification type': 'Unknown', 'All scores': '58.1937302663689;45.8819242017191;36.975033213608', 'All sequences': 'ELMKLAPR;GIFPLGTPR;EKAIFPPR', 'All modified sequences': '_ELMKLAPR_;_GIFPLGTPR_;_EKAIFPPR_', 'id': '6392', 'Protein group IDs': '4453', 'Peptide ID': '6020', 'Mod. peptide ID': '6018', 'Evidence ID': '6357', 'Deamidation (NQ) site IDs': '', 'Oxidation (M) site IDs': ''}, rescoring_features={})

Using the correct spectrum_id_pattern, we can parse all spectra:

[13]:

from ms2rescore.parse_spectra import MSDataType, add_precursor_values, annotate_spectra

_ = add_precursor_values(
    psm_list,
    {MSDataType.ms2_spectra},
    spectrum_path="20161213_NGHF_DBJ_SA_Exp3A_HeLa_1ug_7min_15000_02.mgf",
    spectrum_id_pattern=r".*scan=(\d+)$",
)

annotate_spectra(psm_list, fragmentation_model="cidhcd", ms2_tolerance=20, ms2_tolerance_mode="ppm")

INFO:ms2rescore.parse_spectra:Parsing precursor info from spectrum files...













INFO:ms2rescore.parse_spectra:Annotating MS2 spectra based on search engine identifications...





Now the PSM list contains the annotated MS2 spectrum with the matching spectrum ID, and the spectrum can be accessed as follows:




[14]:





example_psm.spectrum









[14]:







AnnotatedMS2Spectrum(identifier='mzspec=20161213_NGHF_DBJ_SA_Exp3A_HeLa_1ug_7min_15000_02.raw: controllerType=0 controllerNumber=1 scan=2', peaks=224, annotated=13)









Basic features


Basic features from BasicFeatureGenerator can be added from any PSM list and contain simple feature, such as the search engine score, charge state, and absolute precursor mass error.




[15]:





basic_fgen = BasicFeatureGenerator()
basic_fgen.add_features(psm_list)















INFO:ms2rescore.feature_generators.basic:Adding basic features to PSMs.







[16]:





pd.DataFrame(list(psm_list["rescoring_features"]))









[16]:











  
    

      
      charge_n
      charge_1
      charge_2
      charge_3
      charge_4
      charge_5
      charge_6
      abs_ms1_error_ppm
      search_engine_score
      theoretical_mass
      experimental_mass
      mass_error
      pep_len
    

  
  
    

      0
      2
      0
      1
      0
      0
      0
      0
      1.137810
      107.6600
      955.546177
      955.545087
      -0.001090
      11
    

    

      1
      3
      0
      0
      1
      0
      0
      0
      0.598941
      107.7400
      2742.435725
      2742.434081
      -0.001644
      25
    

    

      2
      3
      0
      0
      1
      0
      0
      0
      0.573291
      137.8900
      2870.530688
      2870.529041
      -0.001647
      26
    

    

      3
      2
      0
      1
      0
      0
      0
      0
      0.929236
      22.6410
      1168.559595
      1168.558507
      -0.001088
      14
    

    

      4
      2
      0
      1
      0
      0
      0
      0
      0.666129
      89.4030
      1645.734325
      1645.733227
      -0.001098
      17
    

    

      ...
      ...
      ...
      ...
      ...
      ...
      ...
      ...
      ...
      ...
      ...
      ...
      ...
      ...
    

    

      35551
      2
      0
      1
      0
      0
      0
      0
      0.782088
      76.1000
      1403.720846
      1403.719747
      -0.001099
      12
    

    

      35552
      4
      0
      0
      0
      1
      0
      0
      1.147761
      114.6400
      1906.996328
      1906.994134
      -0.002193
      16
    

    

      35553
      2
      0
      1
      0
      0
      0
      0
      0.813101
      43.1240
      1356.697652
      1356.696547
      -0.001105
      11
    

    

      35554
      3
      0
      0
      1
      0
      0
      0
      0.719097
      23.5620
      2276.122890
      2276.121251
      -0.001639
      18
    

    

      35555
      3
      0
      0
      1
      0
      0
      0
      1.083807
      9.1885
      1513.729834
      1513.728191
      -0.001644
      11
    

  



35556 rows × 13 columns






Now the PSMs contain a simple set of features.


Note that the charge state is present as both an integer and a one-hot encoded vector. This is because charge state can act as both a categorical and a numerical feature. Simply put, sometimes an effect increases or decreases with increasing charge state, and sometimes it is simply different for different charge states.






MS²PIP features


MS²PIP is a machine learning tool that predicts the MS2 spectrum of a peptide given its sequence. It is previously identified MS2 spectra and their corresponding peptide sequences. Because MS²PIP uses the highly performant - but traditional - machine learning approach XGBoost, it can already produce accurate predictions even if trained on smaller spectral libraries. This makes MS²PIP a very flexible platform to train new models on custom datasets. Nevertheless, MS²PIP comes with several
pre-trained models, which we will use in this tutorial.


Because traditional proteomics search engines do not fully consider MS2 peak intensities in their scoring functions, adding rescoring features derived from spectrum prediction tools has proved to be a very effective way to further improve the sensitivity of peptide-spectrum matching. Generating features from MS²PIP predictions follows these steps:


    

  1. Predict MS2 spectra for all PSMs in the dataset, including decoy PSMs and low scoring target PSMs.
    2. 

  2. Annotate the observed MS2 spectra from the original spectrum files (already done before by MS²Rescore)
    3. 

  3. Compare the predicted and observed spectra using various similarity metrics, which are returned as rescoring features.
    4. 





Configuring MS²PIP


MS²PIP requires two parameters:



Note that fragment mass tolerance is no longer configured on MS²PIP directly — it is set once in annotate_spectra (see above) and the resulting AnnotatedMS2Spectrum objects are reused here.




[17]:





ms2pip_fgen = MS2PIPFeatureGenerator(
    model="HCD",
)











Generating MS²PIP features




[18]:





ms2pip_fgen.add_features(psm_list)















INFO:ms2rescore.feature_generators.ms2pip:Adding MS²PIP-derived features to PSMs.
INFO:ms2pip.core:Processing spectra and peptides...
WARNING:ms2pip._utils.psm_input:Skipped 2/35556 PSMs with invalid peptidoforms: 2 unsupported amino acid 'U'







[19]:





pd.DataFrame(list(psm_list["rescoring_features"]))









[19]:











  
    

      
      charge_n
      charge_1
      charge_2
      charge_3
      charge_4
      charge_5
      charge_6
      abs_ms1_error_ppm
      search_engine_score
      theoretical_mass
      ...
      ionb_abs_diff_Q1_norm
      iony_mse
      min_abs_diff_iontype
      abs_diff_Q3
      ionb_abs_diff_Q2_norm
      cos_iony_norm
      iony_std_abs_diff
      dotprod_iony
      abs_diff_Q1_norm
      ionb_max_abs_diff_norm
    

  
  
    

      0
      2
      0
      1
      0
      0
      0
      0
      1.137810
      107.6600
      955.546177
      ...
      0.101514
      0.000048
      0.0
      0.009766
      0.734400
      0.995659
      0.004302
      0.001587
      0.165929
      1.684827
    

    

      1
      3
      0
      0
      1
      0
      0
      0
      0.598941
      107.7400
      2742.435725
      ...
      0.029784
      0.000031
      0.0
      0.003249
      0.264225
      0.998034
      0.004661
      0.005130
      0.012931
      2.034264
    

    

      2
      3
      0
      0
      1
      0
      0
      0
      0.573291
      137.8900
      2870.530688
      ...
      0.000000
      0.000006
      0.0
      0.001948
      0.139024
      0.997942
      0.001813
      0.001095
      0.001910
      2.130072
    

    

      3
      2
      0
      1
      0
      0
      0
      0
      0.929236
      22.6410
      1168.559595
      ...
      0.078716
      0.000448
      0.0
      0.014517
      0.338768
      0.955242
      0.016877
      0.000887
      0.086232
      3.987070
    

    

      4
      2
      0
      1
      0
      0
      0
      0
      0.666129
      89.4030
      1645.734325
      ...
      0.006264
      0.000505
      0.0
      0.010269
      0.162389
      0.982892
      0.019758
      0.000698
      0.031096
      3.818113
    

    

      ...
      ...
      ...
      ...
      ...
      ...
      ...
      ...
      ...
      ...
      ...
      ...
      ...
      ...
      ...
      ...
      ...
      ...
      ...
      ...
      ...
      ...
    

    

      35551
      2
      0
      1
      0
      0
      0
      0
      0.782088
      76.1000
      1403.720846
      ...
      0.072093
      0.000628
      0.0
      0.018626
      0.106503
      0.977440
      0.015875
      0.001575
      0.131833
      1.677736
    

    

      35552
      4
      0
      0
      0
      1
      0
      0
      1.147761
      114.6400
      1906.996328
      ...
      0.001098
      0.000005
      0.0
      0.002591
      0.010763
      0.997821
      0.001660
      0.000832
      0.005419
      3.105403
    

    

      35553
      2
      0
      1
      0
      0
      0
      0
      0.813101
      43.1240
      1356.697652
      ...
      0.096080
      0.000517
      0.0
      0.014355
      0.323174
      0.961556
      0.016832
      0.000125
      0.362248
      3.921681
    

    

      35554
      3
      0
      0
      1
      0
      0
      0
      0.719097
      23.5620
      2276.122890
      ...
      0.000000
      0.000064
      0.0
      0.004335
      0.000000
      0.975086
      0.006132
      0.000309
      0.000000
      3.860645
    

    

      35555
      3
      0
      0
      1
      0
      0
      0
      1.083807
      9.1885
      1513.729834
      ...
      0.910986
      0.000033
      0.0
      0.005959
      1.517434
      0.983124
      0.004654
      0.000018
      0.699540
      2.873089
    

  



35556 rows × 84 columns












Generating DeepLC features




[20]:





deeplc_fgen = DeepLCFeatureGenerator(
    processes=8,
    deeplc_retrain=False,
)
deeplc_fgen.add_features(psm_list)















INFO:ms2rescore.feature_generators.deeplc:Adding DeepLC-derived features to PSMs.
INFO:ms2rescore.feature_generators.deeplc:Predicting retention times with DeepLC...





































INFO:ms2rescore.feature_generators.deeplc:Selecting best head and calibrating predicted retention times per run...









Generating features from the MS2 spectra


The MS2 feature generator scores the pre-annotated spectra against their PSM peptidoforms, computing features such as the fraction of explained intensity, the number and percentage of matched b/y ions, the longest consecutive ion series, and a hyperscore. It reuses the AnnotatedMS2Spectrum objects attached to each PSM by annotate_spectra — no additional annotation is performed here.


The fragmentation_model parameter controls which ion series are scored (e.g. "cidhcd" activates a, b, and y ions).




[21]:





ms2_fgen = MS2FeatureGenerator(
    fragmentation_model="cidhcd",
)
ms2_fgen.add_features(psm_list)















INFO:ms2rescore.feature_generators.ms2:Adding MS2-derived features to PSMs.









Assessing individual feature performance before rescoring




[22]:





from ms2rescore.report.charts import (
    calculate_feature_qvalues,
    feature_ecdf_auc_bar,
    fdr_plot,
    ms2pip_correlation,
    rt_scatter,
    rt_distribution_baseline,
)

features = pd.DataFrame(list(psm_list["rescoring_features"]))







The following plot shows the number of identified PSM at varying FDR thresholds, where the red line indicates the commonly used 1% FDR threshold.




[23]:





psm_list.calculate_qvalues()  # msms.txt does not contain q-values
fig = fdr_plot(psm_list.to_dataframe(), fdr_thresholds=None, log=False)
fig.show()















                                
                    




We can make a similar plot for each of the rescoring features by simply calculating q-values as if the feature were a search engine score.




[24]:





feature_qvalues, feature_ecdf_auc = calculate_feature_qvalues(
    features=features,
    is_decoy=psm_list["is_decoy"],
)









[25]:





feature_names = {
    "basic": basic_fgen.feature_names,
    "ms2pip": ms2pip_fgen.feature_names,
    "deeplc": deeplc_fgen.feature_names,
    "ms2": ms2_fgen.feature_names,
}
feature_generator_map = {name: gen for gen, f_list in feature_names.items() for name in f_list}

feature_qvalues_melt = feature_qvalues.melt(var_name="feature", value_name="q-value")
feature_qvalues_melt["feature_generator"] = feature_qvalues_melt["feature"].map(feature_generator_map)









[26]:





# Plot takes 100MB in saved file...

# fig = px.ecdf(
#     data_frame=feature_qvalues_melt,
#     x="q-value",
#     color="feature_generator",
#     line_dash="feature",
# )
# fig.show()







The larger the area under the curve (AUC) is for each line, the better. We can therefore visualize the same information in a more easy-to-interpret bar chart with the AUCs, which were already calculated by the calculate_feature_qvalues function.




[27]:





feature_ecdf_auc["feature_generator"] = feature_ecdf_auc["feature"].map(feature_generator_map)









[28]:





fig = feature_ecdf_auc_bar(feature_ecdf_auc)
fig.show()















                                
                    




Importantly, a low q-value ECDF AUC does not necessarily mean that a feature is not valuable for rescoring. The more orthogonal a feature is to the other information already used in PSM scoring, the more it will help to confidently identify PSMs. For instance, retention time features are completely independent from the information used in a traditional search engine. When combined with other scoring information, they can provide a substantial boost in sensitivity.


The quality of specific rescoring features can also be assessed in a manner specific to the feature. For instance, we can plot the distribution of the MS²PIP pearson correlations for target PSMs that initially passed the FDR threshold:




[29]:





fig = ms2pip_correlation(
    features=features,
    is_decoy=psm_list["is_decoy"],
    qvalue=psm_list["qvalue"]
)
fig.show()















                                
                    




For DeepLC, we can plot the predicted retention time versus the observed retention time:




[30]:





fig = rt_scatter(
    df=features[(~psm_list["is_decoy"]) & (psm_list["qvalue"] <= 0.01)],  # noqa: E712
    predicted_column="predicted_retention_time_best",
    observed_column="observed_retention_time_best",
)
fig.show()















                                
                    




DeepLC also provides a function to plot the current relative mean absolute error (rMAE) of the predicted retention times against a distribution of 460 benchmark datasets. The lower the rMAE, the better.




[31]:





fig = rt_distribution_baseline(
    df=features[(~psm_list["is_decoy"]) & (psm_list["qvalue"] <= 0.01)],  # noqa: E712
    predicted_column="predicted_retention_time_best",
    observed_column="observed_retention_time_best",
)
fig.show()















                                
                    










Rescoring




Using ristretto directly


MS²Rescore uses ristretto as rescoring engine: a lean, dependency-light (numpy/scikit-learn/pandas only) reimplementation of the Percolator/Käll semi-supervised algorithm, purpose-built to serve MS²Rescore. ristretto takes a plain pandas.DataFrame of identifiers and features, so we first build that DataFrame from the PSM list.




[32]:





import ristretto

feature_names_flat = sorted({f for psm in psm_list for f in psm.rescoring_features})

features_df = pd.DataFrame(
    {
        "spectrum_id": psm_list["spectrum_id"],
        "run": psm_list["run"],
        "is_decoy": psm_list["is_decoy"],
        "peptidoform": [str(p) for p in psm_list["peptidoform"]],
        "score": psm_list["score"],
    }
)
feature_values = (
    pd.DataFrame(list(psm_list["rescoring_features"]))[feature_names_flat]
    .astype(float)
    .fillna(0.0)
)
features_df = pd.concat([features_df, feature_values], axis=1)
features_df.head()









[32]:











  
    

      
      spectrum_id
      run
      is_decoy
      peptidoform
      score
      abs_diff_Q1
      abs_diff_Q1_norm
      abs_diff_Q2
      abs_diff_Q2_norm
      abs_diff_Q3
      ...
      rt_diff_best
      search_engine_score
      spec_mse
      spec_mse_norm
      spec_pearson
      spec_pearson_norm
      spec_spearman
      std_abs_diff
      std_abs_diff_norm
      theoretical_mass
    

  
  
    

      0
      4703
      20161213_NGHF_DBJ_SA_Exp3A_HeLa_1ug_7min_15000_02
      False
      AAAAAAALQAK/2
      107.660
      0.000594
      0.165929
      0.002967
      0.666364
      0.009766
      ...
      0.142506
      107.660
      0.000251
      0.739558
      0.953128
      0.954095
      0.918799
      0.013420
      0.524816
      955.546177
    

    

      1
      13572
      20161213_NGHF_DBJ_SA_Exp3A_HeLa_1ug_7min_15000_02
      False
      [UNIMOD:Acetylation]-AAAAAEQQQFYLLLGNLLSPDNVVR/3
      107.740
      0.000009
      0.012931
      0.000820
      0.319576
      0.003249
      ...
      0.252993
      107.740
      0.000159
      0.520783
      0.937566
      0.959657
      0.882700
      0.011488
      0.527856
      2742.435725
    

    

      2
      13366
      20161213_NGHF_DBJ_SA_Exp3A_HeLa_1ug_7min_15000_02
      False
      [UNIMOD:Acetylation]-AAAAAEQQQFYLLLGNLLSPDNVVRK/3
      137.890
      0.000001
      0.001910
      0.000571
      0.311654
      0.001948
      ...
      0.022973
      137.890
      0.000136
      0.433370
      0.919416
      0.945299
      0.890340
      0.011077
      0.496484
      2870.530688
    

    

      3
      505
      20161213_NGHF_DBJ_SA_Exp3A_HeLa_1ug_7min_15000_02
      False
      AAAAAQGGGGGEPR/2
      22.641
      0.000062
      0.086232
      0.001181
      0.811038
      0.014517
      ...
      1.728728
      22.641
      0.000461
      5.768089
      0.523585
      0.529264
      0.503597
      0.017876
      1.797819
      1168.559595
    

    

      4
      6589
      20161213_NGHF_DBJ_SA_Exp3A_HeLa_1ug_7min_15000_02
      False
      AAAAAWEEPSSGN[UNIMOD:Deamidation]GTAR/2
      89.403
      0.000022
      0.031096
      0.001235
      0.801425
      0.010269
      ...
      0.168237
      89.403
      0.000427
      3.596586
      0.639873
      0.831775
      0.886938
      0.018176
      1.377847
      1645.734325
    

  



5 rows × 124 columns






Then, ristretto.rescore can be called to rescore the PSMs. It runs spectrum-grouped k-fold cross-validation, target-decoy competition, and q-value/PEP calculation, returning a RescoreResult with the PSM-level results plus a peptidoform-level rollup (and peptide-/protein-level rollups, if a peptide_col/protein_col is provided).




[33]:





result = ristretto.rescore(
    features_df,
    run_col="run",
    feature_cols=feature_names_flat,
    model="svm",
    train_fdr=0.01,
)
result.psms.head()















INFO:ristretto.rescore:Rescoring 35556 PSMs (22190 targets, 13366 decoys) with 119 features using model 'svm'
INFO:ristretto.rescore:Fold 1/3: training on 24113 PSMs, scoring 11443
INFO:ristretto.rescore:Fold 2/3: training on 23517 PSMs, scoring 12039
INFO:ristretto.rescore:Fold 3/3: training on 23482 PSMs, scoring 12074
INFO:ristretto.rescore:Spectrum competition: 35556 PSMs -> 10568 best per spectrum
INFO:ristretto.rescore:Computing PSM-level q-values and PEPs
INFO:ristretto.rescore:Rolling up to peptidoform level on column 'peptidoform'







[33]:











  
    

      
      spectrum_id
      run
      is_decoy
      peptidoform
      score
      qvalue
      pep
    

  
  
    

      0
      4703
      20161213_NGHF_DBJ_SA_Exp3A_HeLa_1ug_7min_15000_02
      False
      AAAAAAALQAK/2
      4.254977
      0.000345
      0.00681
    

    

      1
      13572
      20161213_NGHF_DBJ_SA_Exp3A_HeLa_1ug_7min_15000_02
      False
      [UNIMOD:Acetylation]-AAAAAEQQQFYLLLGNLLSPDNVVR/3
      4.000630
      0.000345
      0.00681
    

    

      2
      13366
      20161213_NGHF_DBJ_SA_Exp3A_HeLa_1ug_7min_15000_02
      False
      [UNIMOD:Acetylation]-AAAAAEQQQFYLLLGNLLSPDNVVRK/3
      9.357508
      0.000345
      0.00681
    

    

      4
      6589
      20161213_NGHF_DBJ_SA_Exp3A_HeLa_1ug_7min_15000_02
      False
      AAAAAWEEPSSGN[UNIMOD:Deamidation]GTAR/2
      2.673495
      0.000481
      0.00681
    

    

      12
      1029
      20161213_NGHF_DBJ_SA_Exp3A_HeLa_1ug_7min_15000_02
      False
      AAAGRYR/2
      -1.401264
      0.040268
      0.14212
    

  







We can reuse the fdr_plot function from earlier to plot the number of identifications at each FDR threshold, this time on the rescored PSMs:




[34]:





fig = fdr_plot(result.psms, fdr_thresholds=[0.01])
fig.show()















                                
                    




Or to compare the un-rescored PSMs to the rescored PSMs. ristretto.evaluate runs the exact same competition/q-value/PEP machinery as rescore, but on an already-computed score column instead of training a new model — here, the PSMs’ original search engine score — so the two are directly comparable. We can then use charts.fdr_plot_comparison, the same chart MS²Rescore’s own HTML report uses for this comparison:




[35]:





from ms2rescore.report.charts import fdr_plot_comparison

before_result = ristretto.evaluate(features_df, run_col="run")

fig = fdr_plot_comparison(before_result, result, fdr_threshold=0.01)
fig.show()















INFO:ristretto.rescore:Spectrum competition: 35556 PSMs -> 10568 best per spectrum
INFO:ristretto.rescore:Computing PSM-level q-values and PEPs
INFO:ristretto.rescore:Rolling up to peptidoform level on column 'peptidoform'













                                
                    




Using ristretto directly also allows us to quickly rescore with different feature sets, for instance:


    

  • Only basic features
    • 

  • Basic + MS2 features
    • 

  • Basic + MS2 + MS²PIP features
    • 

  • Basic + MS2 + MS²PIP + DeepLC features
    • 





[36]:





feature_sets = {
    "basic": feature_names["basic"],
    "basic_ms2": feature_names["basic"] + feature_names["ms2"],
    "basic_ms2_ms2pip": feature_names["basic"] + feature_names["ms2"] + feature_names["ms2pip"],
    "basic_ms2_ms2pip_deeplc": feature_names["basic"] + feature_names["ms2"] + feature_names["ms2pip"] + feature_names["deeplc"],
}









[37]:





import plotly.express as px

plot_rows = [
    before_result.psms.sort_values("score", ascending=False)
    .drop_duplicates(subset=["run", "spectrum_id"], keep="first")
    .loc[lambda df: ~df["is_decoy"], ["qvalue"]]
    .assign(feature_set="before rescoring")
]
for name, cols in feature_sets.items():
    result_i = ristretto.rescore(
        features_df, run_col="run", feature_cols=cols, model="svm", train_fdr=0.01
    )
    best = result_i.psms.sort_values("score", ascending=False).drop_duplicates(
        subset=["run", "spectrum_id"], keep="first"
    )
    plot_rows.append(best.loc[~best["is_decoy"], ["qvalue"]].assign(feature_set=name))

fig = px.ecdf(pd.concat(plot_rows), x="qvalue", color="feature_set", log_x=True, ecdfnorm=None)
fig.update_layout(xaxis_title="FDR threshold", yaxis_title="Identified spectra")
fig.show()















INFO:ristretto.rescore:Rescoring 35556 PSMs (22190 targets, 13366 decoys) with 13 features using model 'svm'
INFO:ristretto.rescore:Fold 1/3: training on 24113 PSMs, scoring 11443
INFO:ristretto.rescore:Fold 2/3: training on 23517 PSMs, scoring 12039
INFO:ristretto.rescore:Fold 3/3: training on 23482 PSMs, scoring 12074
INFO:ristretto.rescore:Spectrum competition: 35556 PSMs -> 10568 best per spectrum
INFO:ristretto.rescore:Computing PSM-level q-values and PEPs
INFO:ristretto.rescore:Rolling up to peptidoform level on column 'peptidoform'
INFO:ristretto.rescore:Rescoring 35556 PSMs (22190 targets, 13366 decoys) with 42 features using model 'svm'
INFO:ristretto.rescore:Fold 1/3: training on 24113 PSMs, scoring 11443
INFO:ristretto.rescore:Fold 2/3: training on 23517 PSMs, scoring 12039
INFO:ristretto.rescore:Fold 3/3: training on 23482 PSMs, scoring 12074
INFO:ristretto.rescore:Spectrum competition: 35556 PSMs -> 10568 best per spectrum
INFO:ristretto.rescore:Computing PSM-level q-values and PEPs
INFO:ristretto.rescore:Rolling up to peptidoform level on column 'peptidoform'
INFO:ristretto.rescore:Rescoring 35556 PSMs (22190 targets, 13366 decoys) with 113 features using model 'svm'
INFO:ristretto.rescore:Fold 1/3: training on 24113 PSMs, scoring 11443
INFO:ristretto.rescore:Fold 2/3: training on 23517 PSMs, scoring 12039
INFO:ristretto.rescore:Fold 3/3: training on 23482 PSMs, scoring 12074
INFO:ristretto.rescore:Spectrum competition: 35556 PSMs -> 10568 best per spectrum
INFO:ristretto.rescore:Computing PSM-level q-values and PEPs
INFO:ristretto.rescore:Rolling up to peptidoform level on column 'peptidoform'
INFO:ristretto.rescore:Rescoring 35556 PSMs (22190 targets, 13366 decoys) with 119 features using model 'svm'
INFO:ristretto.rescore:Fold 1/3: training on 24113 PSMs, scoring 11443
INFO:ristretto.rescore:Fold 2/3: training on 23517 PSMs, scoring 12039
INFO:ristretto.rescore:Fold 3/3: training on 23482 PSMs, scoring 12074
INFO:ristretto.rescore:Spectrum competition: 35556 PSMs -> 10568 best per spectrum
INFO:ristretto.rescore:Computing PSM-level q-values and PEPs
INFO:ristretto.rescore:Rolling up to peptidoform level on column 'peptidoform'













                                
                    








Using MS²Rescore’s full pipeline


The steps above show what happens under the hood. For most use cases, the single ms2rescore.rescore() function call handles the entire workflow – parsing PSMs, generating all configured features, rescoring with ristretto, and writing the output files (and, optionally, the HTML report) – from one configuration dictionary (or JSON/TOML file). See the configuration guide for all available options.




[38]:





from ms2rescore import rescore, parse_configurations

configuration = {
    "ms2rescore": {
        "psm_file": "msms.txt",
        "psm_file_type": "msms",
        "spectrum_path": "20161213_NGHF_DBJ_SA_Exp3A_HeLa_1ug_7min_15000_02.mgf",
        "spectrum_id_pattern": r".*scan=(\d+)$",
        "modification_mapping": {
            "gl": "U:Gln->pyro-Glu",
            "ox": "U:Oxidation",
            "ac": "U:Acetylation",
            "de": "U:Deamidation",
        },
        "fixed_modifications": {"U:Carbamidomethyl": ["C"]},
        "feature_generators": {
            "basic": {},
            "ms2": {"fragmentation_model": "cidhcd"},
            "ms2pip": {"model": "HCD"}
        }
    }
}

parsed_config = parse_configurations([configuration])
rescore(parsed_config)















INFO:ms2rescore.parse_psms:Reading PSMs from PSM file (1/1): 'msms.txt'...
WARNING:ms2rescore.parse_psms:Removed 2 PSMs with invalid amino acids.
INFO:ms2rescore.parse_psms:Found 35554 PSMs, of which 37.59% are decoys.
INFO:ristretto.rescore:Spectrum competition: 35554 PSMs -> 10568 best per spectrum
INFO:ristretto.rescore:Computing PSM-level q-values and PEPs
INFO:ristretto.rescore:Rolling up to peptidoform level on column 'peptidoform'
INFO:ristretto.rescore:Spectrum competition: 35554 PSMs -> 10568 best per spectrum
INFO:ristretto.rescore:Computing PSM-level q-values and PEPs
INFO:ristretto.rescore:Rolling up to peptidoform level on column 'peptidoform'
c:\Users\ralfg\git\ms2rescore-2\.venv\Lib\site-packages\psm_utils\psm_list.py:201: FutureWarning: Downcasting object dtype arrays on .fillna, .ffill, .bfill is deprecated and will change in a future version. Call result.infer_objects(copy=False) instead. To opt-in to the future behavior, set `pd.set_option('future.no_silent_downcasting', True)`
  .fillna(0)  # groupby does not play well with None values
INFO:ms2rescore.parse_psms:Removed 4240 PSMs with rank >= 10.
INFO:ristretto.rescore:Spectrum competition: 31314 PSMs -> 10568 best per spectrum
INFO:ristretto.rescore:Computing PSM-level q-values and PEPs
INFO:ristretto.rescore:Rolling up to peptidoform level on column 'peptidoform'
INFO:ristretto.rescore:Spectrum competition: 31314 PSMs -> 10568 best per spectrum
INFO:ristretto.rescore:Computing PSM-level q-values and PEPs
INFO:ristretto.rescore:Rolling up to peptidoform level on column 'peptidoform'
INFO:ristretto.rescore:Spectrum competition: 31314 PSMs -> 10568 best per spectrum
INFO:ristretto.rescore:Computing PSM-level q-values and PEPs
INFO:ristretto.rescore:Rolling up to peptidoform level on column 'peptidoform'
INFO:ristretto.rescore:Rolling up to peptide level on column 'peptide'
INFO:ms2rescore.core:Found 5610 identified PSMs at 1.00% FDR before rescoring.
INFO:ms2rescore.parse_spectra:Parsing precursor info from spectrum files...





































INFO:ms2rescore.parse_spectra:Annotating MS2 spectra based on search engine identifications...
INFO:ms2rescore.feature_generators.basic:Adding basic features to PSMs.
INFO:ms2rescore.feature_generators.ms2:Adding MS2-derived features to PSMs.
INFO:ms2rescore.feature_generators.ms2pip:Adding MS²PIP-derived features to PSMs.
INFO:ms2pip.core:Processing spectra and peptides...
INFO:ms2rescore.core:Rescoring 31314 PSMs with ristretto...
INFO:ristretto.rescore:Spectrum competition: 31314 PSMs -> 10568 best per spectrum
INFO:ristretto.rescore:Computing PSM-level q-values and PEPs
INFO:ristretto.rescore:Rolling up to peptidoform level on column 'peptidoform'
INFO:ristretto.rescore:Spectrum competition: 31314 PSMs -> 10568 best per spectrum
INFO:ristretto.rescore:Computing PSM-level q-values and PEPs
INFO:ristretto.rescore:Rolling up to peptidoform level on column 'peptidoform'
INFO:ristretto.rescore:Rescoring 31314 PSMs (20053 targets, 11261 decoys) with 113 features using model 'svm'
INFO:ristretto.rescore:Fold 1/3: training on 20910 PSMs, scoring 10404
INFO:ristretto.rescore:Fold 2/3: training on 20824 PSMs, scoring 10490
INFO:ristretto.rescore:Fold 3/3: training on 20894 PSMs, scoring 10420
INFO:ristretto.rescore:Computing PSM-level q-values and PEPs
INFO:ristretto.rescore:Rolling up to peptidoform level on column 'peptidoform'
INFO:ristretto.rescore:Rolling up to peptide level on column 'peptide'
INFO:ristretto.rescore:Spectrum competition: 31314 PSMs -> 10568 best per spectrum
INFO:ristretto.rescore:Computing PSM-level q-values and PEPs
INFO:ristretto.rescore:Rolling up to peptidoform level on column 'peptidoform'
INFO:ristretto.rescore:Rolling up to peptide level on column 'peptide'
c:\Users\ralfg\git\ms2rescore-2\.venv\Lib\site-packages\psm_utils\psm_list.py:201: FutureWarning: Downcasting object dtype arrays on .fillna, .ffill, .bfill is deprecated and will change in a future version. Call result.infer_objects(copy=False) instead. To opt-in to the future behavior, set `pd.set_option('future.no_silent_downcasting', True)`
  .fillna(0)  # groupby does not play well with None values
INFO:ms2rescore.core:Identified +1545 (27.54%) PSMs at 1.00% FDR after rescoring, compared to before.
INFO:ms2rescore.core:Writing output to msms.ms2rescore.tsv...