Feature request: built-in converter for Sage (`results.sage.tsv`)

[neely]

Hi MSstatsConvert team,

I'd like to request a built-in cleaning/converter function for Sage, the open-source Rust-based proteomics search engine, similar to the existing support for MaxQuant, OpenMS, Spectronaut, Skyline, etc.

Sage is a fast, open-source DDA/DIA search engine that outputs PSMs as a tab-separated file, typically results.sage.tsv. It's increasingly used as a MaxQuant/MSFragger alternative, so native MSstatsClean/*toMSstatsFormat support would help labs standardizing on Sage move directly into MSstats/MSstatsTMT workflows. github

Representative header and row from results.sage.tsv:

psm_id  peptide  proteins  protein_groups  num_proteins  num_protein_groups  filename  scannr  rank  label  expmass  calcmass  charge  peptide_len  missed_cleavages  semi_enzymatic  isotope_error  precursor_ppm  fragment_ppm  hyperscore  delta_next  delta_best  rt  aligned_rt  predicted_rt  delta_rt_model  ion_mobility  predicted_mobility  delta_mobility  matched_peaks  longest_b  longest_y  longest_y_pct  matched_intensity_pct  scored_candidates  poisson  sage_discriminant_score  posterior_error  spectrum_q  peptide_q  protein_q  protein_group_q  ms2_intensity

Proposed mapping to the standard MSstats long format described in the "Working with MSstatsConvert" vignette:

MSstats column	Sage source
ProteinName	proteins (or protein_groups for grouped analysis)
PeptideSequence	peptide (strip inline mod tags like C[+57.021465] into a separate modification representation, similar to how MaxQuant's Modifications column is handled)
PrecursorCharge	charge
FragmentIon, ProductCharge	not available — fill with NA via columns_to_fill, as already supported
Run	filename
Intensity	ms2_intensity (or a quant column if run alongside a Sage quant step)
IsotopeLabelType	constant "L" for label-free Sage runs

Useful filtering columns already present in Sage output that could plug into the existing score_filtering/pattern_filtering mechanism:

• spectrum_q, peptide_q, protein_q - q-value based filtering (similar to score_filtering with direction = "smaller").
• label - decoy flag (label = -1 for decoys), useful for exact_filtering/pattern_filtering to drop decoy hits before statistical modeling.
• rank - could be used to keep only rank-1 PSMs per spectrum.

Following the existing pattern in the package (built-in MSstatsClean support plus a SagetoMSstatsFormat/SagetoMSstatsTMTFormat wrapper, analogous to MaxQtoMSstatsFormat etc.), it would be great to have:

1. Native MSstatsClean(type = "MSstats", tool = "Sage") support that maps the columns above.
2. Documentation of required/optional Sage columns, mirroring the "Importing and cleaning data" and "Preprocessing" sections of the MSstatsConvert vignette.
3. Guidance on decoy (label) and q-value filtering defaults, since Sage reports both target and decoy PSMs in the same TSV by default.

I can provide a small results.sage.tsv (with a handful of PSMs, including decoys and one or two modified peptides) plus a minimal annotation file if that's helpful for building/testing the converter.

Thanks for considering this, and happy to test a draft implementation,
Ben

results.sage.tsv

[swaraj-neu]

Hi @neely,

Thanks for the detailed proposal, it made this easy to scope. One thing before I start.

In the Sage source, ms2_intensity is the sum of matched b/y ion intensities (scoring.rs: summed_b + summed_y) and feeds the discriminant model, so it reads as a PSM quality feature rather than a quant channel. DOCS.md describes it differently but the code is clear. You hedged on this in your original post and I think you were right to.

So: should the converter read lfq.tsv or results.sage.tsv?

• lfq.tsv has the MS1 areas, which is what MSstats wants. It is wide, one column per file, so it needs a melt. One caveat: with the default combine_charge_states = true, charge is written as -1, so PrecursorCharge would be meaningless unless users turn that off.
• results.sage.tsv is PSM level. In your sample, 25,277 PSMs collapse to 18,358 peptide-charge pairs within a single run. If we go this way, how should duplicates be collapsed: max, sum, or best-scoring PSM?

Defaults I am planning, correct me if any are wrong: filter on spectrum_q at 0.01, keep rank == 1, drop label == -1, ProteinName from proteins, and keep inline modification tags in PeptideSequence. I am avoiding protein_q as a default since shared peptides keep it at 1.0.

One request: a small multi-run example with its annotation file, if you have one. Your sample has a single run, so the annotation join cannot be exercised.

Thanks,
Swaraj

[neely]

Swaraj,

Big thanks for writing so quickly and for moving forward on this. I think your read on scoring.rs is right. Use lfq.tsv, not results.sage.tsv. It has the MS1 areas MSstats actually wants, and ms2_intensity being a discriminant feature rather than a quant channel means results.sage.tsv is the wrong source regardless of collapsing strategy.

On charge: my sample used the default combine_charge_states = true, so charge is -1 throughout, that's why lfq.tsv looked that way. I'll rerun with it set to false and send you an updated file so PrecursorCharge is actually meaningful. Worth handling -1 gracefully in the converter anyway, since plenty of users like me will run with the default on.

Overall, your defaults look right to me (keep in mind, not an expert or the author, but have access to the docs and some user experience): spectrum_q at 0.01, rank == 1, dropping label == -1, ProteinName from proteins, and keeping inline mod tags in PeptideSequence. Skipping protein_q as a default makes sense too, agreed on the shared peptide issue.

I'm putting together a 4v4 run now. I'll send the multi-run lfq.tsv plus an annotation file (Run, Condition, BioReplicate) so you can test the join logic properly.

So that is two files, the one with the combine charge states off, and one as a big 4v4.

Thanks,
Ben

[neely]

Turns out running multiple files with combine_charge_states=false makes things really slow as all charge states are compared.

Here are files from that same single run with charge states. https://we.tl/t-n5T3UV2PAOwTjwFn

Then I did a 4 v 4 run, with an annotation table. I tried doing this without the combined charge states but it was taking too long. So here are those files. https://we.tl/t-F3TW8OqXWyVw7Hww

(I used these links which are good for 3 days, cause are a bit large here)

I guess we will need to be clear about what files to use for MSstats and any required settings.

[swaraj-neu]

Thanks Ben, both sets came through. A few findings from the 4v4 file, two worth your attention.

q_value is not filtered before writing. Values run from 0.0008 to 0.288, and 61% of rows sit above 0.01, even though lfq_settings.peptide_q_value is set to 0.01. That setting seems to control which PSMs build the LFQ traces rather than what gets written out. So the default 0.01 filter will drop about two thirds of your rows. I will call this out clearly in the documentation.

Zeros mean not quantified. 9.3% of intensity cells are zero, with no NAs anywhere. Your single-run file has none at all, which fits: a peptide only shows a zero when it was quantified in some runs but not others. So I am converting zeros to NA before modelling. Say if you disagree.

Good news: peptide plus charge is unique in both your files, so no aggregation is needed. And there are no decoys in lfq.tsv, so no decoy filter either.

Two small things. Your annotation has an IsotopeLabelType column, which is not needed since LFQ is label-free and the converter sets it to "L" itself. I will drop it and log a message, so your file still works as is. And I am passing charge = -1 through unchanged rather than converting it to NA, with a message explaining that charge states were combined.

Thanks,
Swaraj

[neely]

Swaraj,

Thanks again for quick responses and work, and for pointing those things out. The zero-to-NA conversion is correct and desired, and one request is could you make the Sage related documentation explicitly state that peptide_q_value filtering must be applied by the user post-hoc, since that's the part most likely to trip up other users who assume the file is already FDR-filtered. It seems peptide_q_value in Sage's LFQ settings is a threshold used internally, Sage intentionally writes all quantified peptide-charge rows regardless of their q-value so that downstream tools can apply their own FDR filtering choice rather than having Sage silently drop rows before you ever see them (or that is how it appears by looking at that lfq.tsv). Also, dropping the redundant IsotopeLabelType column is reasonable since it's a TMT/isobaric-labeling artifact column that has no meaning for LFQ data. And last thing, passing charge = -1 through unchanged for combined-charge-state rows is also sensible bookkeeping, since your config has lfq_settings.combine_charge_states = true, so -1 is Sage's sentinel value meaning "this trace merges multiple charge states," and preserving it with a logged explanation is more transparent than silently converting it to NA.

Thanks again,
Ben

[swaraj-neu]

Thanks Ben, adding that section now.

One clarification worth being explicit about: the converter applies the filter for you rather than leaving it entirely post-hoc. filter_with_Qvalue defaults to TRUE with qvalue_cutoff at 0.01, so a default call returns FDR-filtered output. The documentation will state that lfq.tsv arrives unfiltered, that lfq_settings.peptide_q_value is an internal threshold rather than an output filter, what the defaults do, and the 23 to 61 percent row reduction observed on your two files so the scale is clear.

Good to have the other three confirmed, and all are already implemented as described.

[tonywu1999]

@neely

Swaraj is one vacation for the next 2 weeks, so I'll be around to follow up on some additional questions.

One question I wanted to confirm, it might've been answered already in the earlier comments: Is there any report from Sage for DIA that reports an MS2 intensity per fragment?

Typically when working with DIA data, MSstats uses the MS2 intensity per fragment as the quantifications of interest (DIANN and Spectronaut both report this). My intuition is that MS1 is susceptible to a lot of interference from other peptides, which is why MS2 is preferred (there are some edge cases MS1 may be better, I'll make that out of scope for now). Meanwhile, a sum of MS2 intensities is not sufficient because sums are not robust to outliers (imagine summing 1, 1, 1, 1, and 5000 as an outlier, which leads to a sum of 5004).

My gut is telling me we can easily set up a converter for DDA in Sage (as we use MS1 here), but DIA seems tricky right now.

Tony