results.sage.pin is missing flanking amino acids and produces an error in Percolator with the -f flag

[jesseuwu]

We are testing out Sage in some of our experiments and are trying to pair it with Percolator using the -f flag to output proteins.

After running Sage, we run Percolator with the command:

percolator -f db.fasta -P "reverse_" -l out.file results.sage.pin

We are getting this error in Percolator:

ERROR: Reading tab file, the peptide sequence XXXXXXXXXX with PSM id 12345 does not contain one or two of its flanking amino acids.

From looking at the results.sage.pin file, the Peptide column contains XXXXXXXXXX but not the flanking amino acids, e.g., K.XXXXXXXXXX.R.

[lazear]

Idk if Sage has ever been tested against percolator (I could never get it run without lots of segfaults). Maybe try mokapot? Or just fudging and adding "K" and "R" and making sure it completes?

In general, Sage's LDA is pretty competitive against mokapot and percolator (within -5 to +5% PSMs), and the .pin files generated are really a legacy artifact

[jesseuwu]

I was never able to get Mokapot to work myself, it would always throw exceptions no matter what I tried. 😅

Regarding adding amino acids, it was just an example, but I believe the flanking amino acids are actually supposed to be the 2 amino acids around the peptide returned in results.sage.pin, eg, in one of our examples Sage returned the sequence GVVDSDDLPLNVSR from this protein: https://rest.uniprot.org/uniprotkb/P14625.fasta and with the flanking amino acids I think it's supposed to be K.GVVDSDDLPLNVSR.E. I assume adding 2 random amino acids would produce wrong results.

The built-in LDA does seem fairly good but my lab member wanted to compare the results to the tools they've previously used so I figured I'd ask, but I will pass this note onto them.

[lazear]

Yeah I mean I think a python script would suffice to remap sequence -> uniprot and find the flanking AAs. My pushback on the several times this request has come up (or related, like recording sequence positions) is just from the added complexity of handling multiple protein mapping/protein grouping parsimony. You could imagine cases where the same GVVDSDDLPLNVSR peptide maps to multiple proteins (or even the same protein at different locations within its sequence) with different amino acids on the flanks, like GVVDSDDLPLNVSR.E and GVVDSDDLPLNVSR.D.

Since this is something that is super trivial for the few end users who need this functionality to do on their end (and they can make the decision about how to handle multiple proteins) and relatively complex to integrate into Sage, I don't have any plans to directly add this. If you need help writing a python script to do this, I would be happy to help.