Lesson 4.7
Mass Spectrometry Meets Golf
Watch the lectureAdapting Sequencing for Spectra with Errors
From theoretical to real spectra
Although CyclopeptideSequencing successfully reconstructed Tyrocidine B1, this algorithm only works in the case of an ideal spectrum, i.e., when the experimental spectrum of a peptide coincides exactly with its theoretical spectrum. This inflexibility of CyclopeptideSequencing presents a practical barrier, since mass spectrometers generate "noisy" spectra that are far from ideal — they are characterized by having both false masses and missing masses. A false mass is present in the experimental spectrum but absent from the theoretical spectrum; a missing mass is present in the theoretical spectrum but absent from the experimental spectrum.
For example, compare the following theoretical and (simulated) experimental spectra of the cyclic peptide NQEL. Masses that are missing from the experimental spectrum are shown in blue, and false masses in the experimental spectrum are shown in green.
What is particularly worrisome about this example is that the mass of the amino acid E (129) is missing, and the mass of the amino acid V (99) is false; as a result, the first step of CyclopeptideSequencing would establish {V, L, N, Q} as the amino acid composition of our candidate peptides, which is incorrect. In fact, any false or missing mass will cause CyclopeptideSequencing to throw out the correct peptide, because its theoretical spectrum differs from the experimental spectrum.
STOP and Think How would you reformulate the Cyclopeptide Sequencing Problem to handle experimental spectra with errors?
Adapting cyclopeptide sequencing for spectra with errors
To generalize the Cyclopeptide Sequencing Problem to handle noisy spectra, we need to relax the requirement that a candidate peptide’s theoretical spectrum must match the experimental spectrum exactly, and instead incorporate a scoring function that will select the peptide whose theoretical spectrum matches the given experimental spectrum the most closely. Given a cyclic peptide Peptide and a spectrum Spectrum, we define Score(Peptide, Spectrum) as the number of masses shared between Cyclospectrum(Peptide) and Spectrum. Recall our previous example:
If Spectrum is the experimental spectrum in the bottom row of this figure, then Score(NQEL, Spectrum) = 11 (the number of columns shared by the theoretical and experimental spectra, shown in black).
The scoring function should take into account the multiplicities of shared masses, i.e., how many times they occur in each spectrum. For example, suppose that Spectrum is the theoretical spectrum of NQEL; for this spectrum, mass 242 has multiplicity 2. If 242 has multiplicity 1 in the experimental spectrum of Peptide, then 242 contributes 1 to Score(Peptide, Spectrum). If 242 has multiplicity 2 or more in the experimental spectrum of Peptide, then 242 contributes 2 to Score(Peptide, Spectrum)
Cyclopeptide Scoring Problem: Compute the score of a cyclic peptide against a spectrum.
- Input: An amino acid string Peptide and a collection of integers Spectrum.
- Output: The score of Peptide against Spectrum, Score(Peptide, Spectrum).
Code Challenge Solve the Cyclopeptide Scoring Problem at Cogniterra or Rosalind.Solve it with the autograderFree to start in the interactive text
We now can redefine the Cyclopeptide Sequencing Problem for noisy spectra.
Cyclopeptide Sequencing Problem (for spectra with errors): Find a cyclic peptide having maximum score against an experimental spectrum.
- Input: A collection of integers Spectrum.
- Output: A cyclic peptide Peptide maximizing Score(Peptide, Spectrum) over all peptides Peptide with mass equal to ParentMass(Spectrum).
Our goal is to adapt the CyclopeptideSequencing algorithm to find a peptide with maximum score. Remember that this algorithm had a stringent bounding step in which all candidate linear peptides having inconsistent spectra were thrown out. For example, we saw that the linear peptide VKF is inconsistent with the theoretical spectrum of the cyclic peptide Tyrocidine B1. However, we perhaps should not ban VKF in the case of experimental spectra, since they can have missing masses. Thus, we need to revise this step to include more candidate linear peptides, while still ensuring that the number of peptides that we consider does not grow out of control.
STOP and Think How can we limit the growth of the list of candidate linear peptides in the case of experimental spectra?
To limit the number of candidate linear peptides under consideration, we will replace Peptides with Leaderboard, which holds the N highest scoring candidate peptides for further extension. At each step, we will expand all candidate peptides found in Leaderboard, then eliminate those peptides whose newly calculated scores are not high enough to keep them on the Leaderboard. This idea is similar to the notion of a “cut” in a golf tournament; after the cut, only the top N golfers are allowed to play in the next round, since they are the only players who have a reasonable chance of winning.
To be fair, a cut should include anyone who is tied with the Nth-place competitor. Thus, Leaderboard should be trimmed down to the “N highest-scoring linear peptides including ties”, which may include more than N peptides. Given a list of peptides Leaderboard, a spectrum Spectrum, and an integer N, define Trim(Leaderboard, Spectrum, N) as the collection of the top N highest-scoring linear peptides in Leaderboard (including ties) with respect to Spectrum.
STOP and Think Our peptide scoring function currently assumes that peptides are circular, but peptides should technically be scored as linear peptides until the final step. How would you score a linear peptide against a spectrum?
Note that Score(Peptide, Spectrum) currently only scores Peptide against Spectrum if Peptide is cyclic. However, to generalize this scoring function when Peptide is linear, we simply exclude those subpeptides of Peptide that wrap around the end of the string, resulting in a function LinearScore(Peptide, Spectrum). For example, if Spectrum is the experimental spectrum of NQEL, then you can verify that LinearScore(NQEL, Spectrum) = 8.
We now introduce LeaderboardCyclopeptideSequencing.
LeaderboardCyclopeptideSequencing(Spectrum, N)
Leaderboard ← set containing only the empty peptide
LeaderPeptide ← empty peptide
while Leaderboard is non-empty
Leaderboard ← Expand(Leaderboard)
for each Peptide in Leaderboard
if Mass(Peptide) = ParentMass(Spectrum)
if Score(Peptide, Spectrum) > Score(LeaderPeptide, Spectrum)
LeaderPeptide ← Peptide
else if Mass(Peptide) > ParentMass(Spectrum)
remove Peptide from Leaderboard
Leaderboard ← Trim(Leaderboard, Spectrum, N)
output LeaderPeptideThe tricky aspect of implementing LeaderboardCyclopeptideSequencing is making sure that the Trim function is implemented correctly. Check out "Charging Station: Trimming the Peptide Leaderboard" in the print companion or at Cogniterra for some help with trimming Leaderboard.
Code Challenge Implement LeaderboardCyclopeptideSequencing at Cogniterra or at Rosalind.Solve it with the autograderFree to start in the interactive text
- Input: An integer N and a collection of integers Spectrum.
- Output: LeaderPeptide after running LeaderboardCyclopeptideSequencing(Spectrum, N).
We point out that because a linear peptide giving rise to the highest-scoring cyclic peptide may be trimmed from the leaderboard early on, LeaderboardCyclopeptideSequencing is a heuristic, not guaranteed to correctly solve the Cyclopeptide Sequencing Problem. When we develop a heuristic, we must ask ourselves: how accurate is it? Consider the simulated spectrum Spectrum10 of Tyrocidine B1 shown below, which has approximately 10% missing/false masses. Note that the blue masses are not actually in the spectrum, but we show them so that it is clear which masses are missing.
Applying LeaderboardCyclopeptideSequencing to this spectrum (with N = 1000) results in the correct cyclic peptide VKLFPWFNQY, which has a score of 86.
So far, LeaderboardCyclopeptideSequencing has worked well, but as the number of errors increases, so does the likelihood that this algorithm will return an incorrect peptide. Let’s see how this algorithm performs on a noisier simulated spectrum; below, we show Spectrum25 for Tyrocidine B1, which has 25% missing/false masses.
Exercise Break Run LeaderboardCyclopeptideSequencing on Spectrum25 with N = 1000. (Solve this at Cogniterra). You should find 38 linear peptides of maximum score 83 (corresponding to 15 different cyclic peptides). What are they?
When run on Spectrum25, LeaderboardCyclopeptideSequencing (with N = 1000) identifies VKLFPADFNQY (score: 83) as a highest-scoring cyclic peptide instead of the correct peptide VKLFPWFNQY (score: 82). These two peptides are similar, owing to the fact that the combined mass of A (71) and D (115) is equal to the mass of W (186).
STOP and Think How could we have eliminated the incorrect peptide VKLFPADFNQY from consideration for Spectrum25?
Notice that although the correct and incorrect peptides are similar, their amino acid compositions differ. If we could figure out the amino acid composition of Tyrocidine B1 from its spectrum alone and run LeaderboardCyclopeptideSequencing on this smaller alphabet (rather than on the alphabet of all amino acids), then we could eliminate the incorrect peptide VKLFPADFNQY from consideration.