Find a cyclic peptide with a given non-ideal spectrum

Question

For a course I need to build a peptide from scratch that matches the given spectrum. So I start with the amino acids that are in the given spectrum and filter those out that do not have the right criteria. Then I trim this so called leaderboard, and expand them.

Here I add all the amino-acids (20) to the already existing peptides so e.g. a leaderboard of length 4 becomes length 80. After a while I get with a leaderboard that is around 4000. And of these 4000 we again check all the criteria and trim and so on.

I don't think the function CycloPepSeq is the problem but the others might be. Trim trims the Leaderboard, which is built up in CycloPepSeq based on a score that is calculated in Score. Only the first N (with ties) can be used in the next step. Expand obviously expands all the peptides with all the amino-acids available, so I think a for in a for loop is the best I can get? Cyclic_spectrum and linear_spectrum build an ideal spectrum for each peptide.

My brute force code isn't efficient enough to pass the time limit set by the site where I need to hand it in. I get a time limit exceeded for this input: CycloPepSeq(325, (0, 71, 71, 71, 87, 97, 97, 99, 101, 103, 113, 113, 114, 115, 128, 128, 129, 137, 147, 163, 163, 170, 184, 184, 186, 186, 190, 211, 215, 226, 226, 229, 231, 238, 241, 244, 246, 257, 257, 276, 277, 278, 299, 300, 312, 316, 317, 318, 318, 323, 328, 340, 343, 344, 347, 349, 356, 366, 370, 373, 374, 391, 401, 414, 414, 415, 419, 427, 427, 431, 437, 441, 446, 453, 462, 462, 462, 470, 472, 502, 503, 503, 511, 515, 529, 530, 533, 533, 540, 543, 547, 556, 559, 569, 574, 575, 584, 590, 600, 600, 604, 612, 616, 617, 630, 640, 640, 643, 646, 648, 660, 671, 683, 684, 687, 693, 703, 703, 719, 719, 719, 729, 730, 731, 737, 740, 741, 745, 747, 754, 774, 780, 784, 790, 797, 800, 806, 818, 826, 827, 832, 833, 838, 846, 846, 847, 850, 868, 869, 877, 884, 889, 893, 897, 903, 908, 913, 917, 930, 940, 947, 956, 960, 960, 961, 964, 965, 966, 983, 983, 985, 1002, 1009, 1010, 1011, 1021, 1031, 1031, 1036, 1053, 1054, 1058, 1059, 1062, 1063, 1074, 1076, 1084, 1092, 1103, 1113, 1122, 1124, 1130, 1133, 1134, 1145, 1146, 1146, 1149, 1150, 1155, 1156, 1171, 1173, 1174, 1187, 1191, 1193, 1200, 1212, 1221, 1233, 1240, 1242, 1246, 1259, 1260, 1262, 1277, 1278, 1283, 1284, 1287, 1287, 1288, 1299, 1300, 1303, 1309, 1311, 1320, 1330, 1341, 1349, 1357, 1359, 1370, 1371, 1374, 1375, 1379, 1380, 1397, 1402, 1402, 1412, 1422, 1423, 1424, 1431, 1448, 1450, 1450, 1467, 1468, 1469, 1472, 1473, 1473, 1477, 1486, 1493, 1503, 1516, 1520, 1525, 1530, 1536, 1540, 1544, 1549, 1556, 1564, 1565, 1583, 1586, 1587, 1587, 1595, 1600, 1601, 1606, 1607, 1615, 1627, 1633, 1636, 1643, 1649, 1653, 1659, 1679, 1686, 1688, 1692, 1693, 1696, 1702, 1703, 1704, 1714, 1714, 1714, 1730, 1730, 1740, 1746, 1749, 1750, 1762, 1773, 1785, 1787, 1790, 1793, 1793, 1803, 1816, 1817, 1821, 1829, 1833, 1833, 1843, 1849, 1858, 1859, 1864, 1877, 1886, 1890, 1893, 1900, 1900, 1903, 1904, 1918, 1922, 1930, 1930, 1931, 1961, 1963, 1971, 1971, 1971, 1980, 1987, 1992, 1996, 2002, 2006, 2006, 2014, 2018, 2019, 2019, 2032, 2042, 2059, 2060, 2063, 2067, 2077, 2084, 2086, 2089, 2090, 2093, 2105, 2110, 2115, 2115, 2116, 2117, 2121, 2133, 2134, 2155, 2156, 2157, 2176, 2176, 2187, 2189, 2192, 2195, 2202, 2204, 2207, 2207, 2218, 2222, 2243, 2247, 2247, 2249, 2249, 2263, 2270, 2270, 2286, 2296, 2304, 2305, 2305, 2318, 2319, 2320, 2320, 2330, 2332, 2334, 2336, 2336, 2346, 2362, 2362, 2362, 2433))

from itertools import islice, cycle
def CycloPepSeq(N, spectrum):
    """
    >>> CycloPepSeq(10, (0, 71, 113, 129, 147, 200, 218, 260, 313, 331, 347, 389, 460))
    (71, 147, 113, 129)
    >>> CycloPepSeq(376, [0, 97, 101, 115, 128, 128, 129, 147, 198, 216, 225, 256, 257, 262, 276, 313, 326, 353, 363, 385, 391, 404, 441, 454, 460, 482, 492, 519, 532, 569, 583, 588, 589, 620, 629, 647, 698, 716, 717, 717, 730, 744, 748, 845, 507, 687])
    (147, 129, 128, 128, 97, 101, 115)
    >>> CycloPepSeq(216, [0, 87, 113, 114, 114, 156, 200, 227, 228, 243, 270, 314, 341, 356, 357, 384, 428, 470, 470, 471, 497, 584, 488])
    (156, 114, 114, 113, 87)
    """
    AAmass = [57, 71, 87, 97, 99, 101, 103, 113, 114, 115, 128, 129, 131, 137, 147, 156, 163, 186]
    LeaderBoard = [(AA, ) for AA in AAmass if AA in spectrum]
    LeaderBoard = expand(LeaderBoard)
    LeaderPeptide = tuple()
    ParentMass = max(spectrum)
    while LeaderBoard:
        partLeaderBoard = list(LeaderBoard)
        for peptide in LeaderBoard:
            pepMass = sum(peptide)
            if pepMass == ParentMass:
                if score(peptide, spectrum, True) > score(LeaderPeptide, spectrum, True):
                    LeaderPeptide = tuple(peptide)
            elif pepMass > ParentMass:
                partLeaderBoard.remove(peptide)
        LeaderBoard = trim(N, spectrum, partLeaderBoard)
        LeaderBoard = expand(LeaderBoard)
    return LeaderPeptide


def trim(N, spectrum, Leaderboard):
    if len(Leaderboard) <= N:
        return Leaderboard
    scores = []
    for pep in Leaderboard:
        scores.append((score(pep, spectrum), pep))
    scores = sorted(scores, reverse=True)
    for j in range(N, len(scores)):
        if scores[j][0] < scores[N-1][0]:
            result = [x[1] for x in scores[:j]]
            return result
    return [x[1] for x in scores]


def expand(LeaderBoard):
    AAmass = [(57,), (71,), (87,), (97,), (99,), (101,), (103,), (113,), (114,), (115,), (128,), (129,), (131,), (137,), (147,), (156,), (163,), (186,)]
    partLeaderBoard = []
    for x in LeaderBoard:
        for y in AAmass:
            partLeaderBoard.append((x + y))
    return partLeaderBoard


def score(peptide, spectrum, cyclic=False):
    spectrum = list(spectrum)
    if cyclic:
        idealspectrum = list(cyclic_spectrum(peptide))
    else:
        idealspectrum = list(linear_spectrum(peptide))
    score = 0
    for i in idealspectrum:
        if i in spectrum:
            score += 1
            spectrum.remove(i)
    return score

def cyclic_spectrum(peptide):
    result = [0]
    for num in range(1, len(peptide)):
        for start in range(len(peptide)):
            group = islice(cycle(peptide), start, start + num)
            result.append(sum(group))
    result.append(sum(peptide))
    return tuple(sorted(result))


def linear_spectrum(peptide):
    PrefixMass = [0]
    lnPep = len(peptide)
    for i in range(lnPep):
        PrefixMass.append(PrefixMass[i] + peptide[i])
    LinearSpectrum = [0]
    for i in range(lnPep):
        for j in range(i + 1, lnPep + 1):
            LinearSpectrum.append(PrefixMass[j] - PrefixMass[i])
    return tuple(sorted(LinearSpectrum))

Answer 1

In the absence of further information, it is difficult to offer advice that addresses the core of the algorithm. In particular, I would want more information about how it's meant to work (and at least half a dozen diverse test cases) before looking for ways to make it scale better with different sizes of inputs.

---

In the mean time, enter cProfile.

I have a main function running CycloPepSeq(10, (0, 71, 113, 129, 147, 200, 218, 260, 313, 331, 347, 389, 460)) 1000 times. Then running

python -m cProfile -s cumtime  peptides_mod.py

Gives a chart of how much time is spent in each function. tottime is the time spent just in the body of the function, cumtime is the time spent in that function and functions it calls.

(147, 71, 129, 113)
         14919007 function calls in 6.673 seconds

   Ordered by: cumulative time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.000    0.000    6.673    6.673 {built-in method builtins.exec}
        1    0.000    0.000    6.673    6.673 peptides_mod.py:1(<module>)
        1    0.001    0.001    6.673    6.673 peptides_mod.py:86(main)
     1000    0.239    0.000    6.672    0.007 peptides_mod.py:2(CycloPepSeq)
   574000    1.335    0.000    5.128    0.000 peptides_mod.py:51(score)
     4000    0.222    0.000    4.834    0.001 peptides_mod.py:26(trim)
   530000    2.192    0.000    2.852    0.000 peptides_mod.py:74(linear_spectrum)
  3425000    0.851    0.000    0.851    0.000 {method 'remove' of 'list' objects}
    44000    0.372    0.000    0.585    0.000 peptides_mod.py:64(cyclic_spectrum)
   577000    0.574    0.000    0.574    0.000 {built-in method builtins.sorted}
  7548000    0.451    0.000    0.451    0.000 {method 'append' of 'list' objects}
  1496000    0.256    0.000    0.256    0.000 {built-in method builtins.sum}
     5000    0.122    0.000    0.175    0.000 peptides_mod.py:40(expand)
   710000    0.050    0.000    0.050    0.000 {built-in method builtins.len}
     1000    0.004    0.000    0.004    0.000 peptides_mod.py:8(<listcomp>)
     3000    0.003    0.000    0.003    0.000 peptides_mod.py:35(<listcomp>)
     1000    0.001    0.000    0.001    0.000 {built-in method builtins.max}
        1    0.000    0.000    0.000    0.000 {built-in method builtins.print}
        1    0.000    0.000    0.000    0.000 <frozen importlib._bootstrap>:989(_handle_fromlist)
        1    0.000    0.000    0.000    0.000 {built-in method builtins.hasattr}
        1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}

So, the whole thing took 6.673 seconds. Of that, 1.335 was spent in score (tottime, so not subfunctions) and 2.192 was spent in linear_spectrum. The built in remove also sits at 0.851, just because removing from an array is quite slow, and generally to be avoided.

Those are the places to look for most obvious improvements.

score

Isolating the following code into its own function reveals that it takes 0.939 seconds tottime, 1.281 cumtime. The difference is entirely because removing an element from a list is quite expensive.

for i in idealspectrum:
    if i in spectrum:
        score += 1
        spectrum.remove(i)

So, looking at what this is doing, it's essentially returning the size of the set intersection between spectrum and idealspectrum. That's easy enough to do. If idealspectrum and spectrum are sets instead of lists it can be written len(idealspectrum & spectrum)

Now score is

def score(peptide, spectrum, cyclic=False):
    spectrum = set(spectrum)
    if cyclic:
        idealspectrum = set(cyclic_spectrum(peptide))
    else:
        idealspectrum = set(linear_spectrum(peptide))

    # Score is the cardinality of the set intersection between the two sets
    return len(idealspectrum & spectrum)

It takes 0.880s, and the whole program is down to 5.8. As a sanity check (although by no means a full test suite) the result has not changed.

As mentioned in the comments, the code in question did not do a set intersection as I had understood, but rather does a multiset intersection. The python object that best models a multiset is collections.Counter. I'll put further measurements at the end.

linear_spectrum

Now, an interesting thing about linear_spectrum is that it is only called by score. So one thing we can do is drop tuple and sorted because the result is going to be put immediately into a set. That shaves another half a second.

We can also replace the creation of PrefixMass with a call to itertools.accumulate. There isn't a faster algorithm for what you're doing, but the built in functions do just tend to run faster, shaving another 0.2 seconds.

Unfortunately the following nested loop is what takes a lot of the time

LinearSpectrum = [0]
for i in range(lnPep):
    for j in range(i + 1, lnPep + 1):
        LinearSpectrum.append(PrefixMass[j] - PrefixMass[i])
return LinearSpectrum

We can shave off a little bit of this, by writing it as list comprehension.

def linear_spectrum(peptide):
    PrefixMass = [0] + list(accumulate(peptide))
    lnPep = len(peptide)

    return [0] + [PrefixMass[pm2] - PrefixMass[pm1] for pm1 in range(lnPep) for pm2 in range(pm1 + 1, lnPep + 1)]

Unfortunately the task for this little function seems to be inherently quadratic, so we probably can't speed it up all that much more without more extensive rearchitecture.

Now the overall function takes 4.746 seconds, which is an noticeable improvement.

Duplicate work

It may be possible to optimise other functions in a similar way, again guided in terms of importance by how much they take up in the profiler.

Of course, the most efficient work is always the work you don't do. For example, we already discovered that we're (still) taking a considerable amount of time inside score. Now, in the topmost function we find this line:

if score(peptide, spectrum, True) > score(LeaderPeptide, spectrum, True):

That recalculates the Leader's score every time there's a contender. Instead, we can do something like this.

    contender_score = score(peptide, spectrum_set, True)
    if contender_score > LeaderPeptideScore:
        LeaderPeptide = tuple(peptide)
        LeaderPeptideScore = contender_score

(LeaderPeptideScore needs initializing once at the start)

The implication here is that every peptide you test, you only test once and keep the result. And of course less work implies less time doing it.

Another point on the less work thing: until now I have been converting spectrum from a tuple to a set in each call to score. Score is run 553 times per call to CycloPepSeq Instead, I can do it once in the root of CycloPepSeq.

Now running at around 4.05 seconds.

Remove

Oh, I forgot that remove is slowing things down. One of the remove calls got diked out in the switch to using set intersection, but there's still one in partLeaderBoard.

It's more efficient to filter things out in one go, and especially while doing something else, than to take them out separately, and especially separately one element at a time.

Here's a way to do that

    augmentedLeaderBoard = [ (p, sp) for (p, sp) in ((peptide, sum(peptide)) for peptide in LeaderBoard) if sp <= ParentMass ]
    partLeaderBoard = [ p for (p, sp) in augmentedLeaderBoard ]
    for peptide, pepMass in augmentedLeaderBoard:

---

Another look at the profile information.

python -m cProfile -s cumtime  peptides_mod.py
(147, 71, 129, 113)
         6780007 function calls in 3.707 seconds

   Ordered by: cumulative time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.000    0.000    3.707    3.707 {built-in method builtins.exec}
        1    0.000    0.000    3.707    3.707 peptides_mod.py:1(<module>)
        1    0.001    0.001    3.707    3.707 peptides_mod.py:83(main)
     1000    0.067    0.000    3.706    0.004 peptides_mod.py:2(CycloPepSeq)
     4000    0.231    0.000    2.740    0.001 peptides_mod.py:32(trim)
   553000    0.626    0.000    2.558    0.000 peptides_mod.py:56(score)
   530000    0.693    0.000    1.616    0.000 peptides_mod.py:76(linear_spectrum)
   530000    0.891    0.000    0.891    0.000 peptides_mod.py:80(<listcomp>)
     4000    0.130    0.000    0.397    0.000 peptides_mod.py:18(<listcomp>)
    23000    0.179    0.000    0.282    0.000 peptides_mod.py:66(cyclic_spectrum)
   940000    0.141    0.000    0.267    0.000 peptides_mod.py:18(<genexpr>)
    26000    0.247    0.000    0.247    0.000 {built-in method builtins.sorted}
  1223000    0.183    0.000    0.183    0.000 {built-in method builtins.sum}
     5000    0.116    0.000    0.165    0.000 peptides_mod.py:46(expand)
  1753000    0.104    0.000    0.104    0.000 {method 'append' of 'list' objects}
  1179000    0.072    0.000    0.072    0.000 {built-in method builtins.len}
     4000    0.019    0.000    0.019    0.000 peptides_mod.py:19(<listcomp>)
     3000    0.003    0.000    0.003    0.000 peptides_mod.py:41(<listcomp>)
     1000    0.001    0.000    0.001    0.000 peptides_mod.py:12(<listcomp>)
     1000    0.001    0.000    0.001    0.000 {built-in method builtins.max}
        1    0.000    0.000    0.000    0.000 {built-in method builtins.print}
        1    0.000    0.000    0.000    0.000 <frozen importlib._bootstrap>:989(_handle_fromlist)
        1    0.000    0.000    0.000    0.000 {built-in method builtins.hasattr}
        1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}

The big thing is that the headline number is down from 6.673 to 3.707 seconds.

There are some new entries; some chunky list comprehensions for example are pulling a bit more weight. Trim remains interesting. It is lightweight in itself, but in terms of cumulative time it is now over 70% of the runtime. Most of that lives in score. (In fact 96% of calls to score are from Trim in the sample input). The best chance at further gains would probably come from thinking about score and Trim together.

---

This is a useful basket of tricks, and top of the basket of course is use a profiler. I'm sure you can apply them further, and perhaps shave the runtime down by another second or two. (Remember to profile your initial code on your own machine: your timings won't match mine)

However, I have a suspicion that this is not going to be enough for you. In general these efficiency tasks aren't won or lost based on a factor of two or three. Instead, they're won or lost on algorithm design, and making sure that you don't unnecessarily duplicate work.

Score gets called 553 times with the input in the code. The last 10 values of peptide it considers are:

(71, 113, 147, 87), 
(71, 113, 147, 97), 
(71, 113, 147, 99), 
(71, 113, 147, 101), 
(71, 113, 147, 103), 
(71, 113, 147, 113), 
(71, 113, 147, 114), 
(71, 113, 147, 115), 
(71, 113, 147, 128), 
(71, 113, 147, 129)

That, to me, screams "Duplicate work". Microoptimsations and library calls like the above will get you some way, but you probably need to look for ways to use the fact that you've already seen these prefixes before.

----------

Completely separate to the above, as less of a performance focussed code review, I noticed the following.

AAmass = [57, 71, 87, 97, 99, 101, 103, 113, 114, 115, 128, 129, 131, 137, 147, 156, 163, 186]

This list is defined twice, once in CycloPepSeq and once in expand. It would be preferable to avoid this, not least to help mitigate the risk that they get out of sync. (Granted one is values and the other is single value tuples, but it's the underlying data that's being unnecessarily repeated.)

---

"""
>>> leaderboard_cyclopeptide_sequencing(10, (0, 71, 113, 129, 147, 200, 218, 260, 313, 331, 347, 389, 460))
(71, 147, 113, 129)
"""

I note that this function name leaderboard_cyclopeptide_sequencing has appeared out of nowhere. I'm guessing it's an old name for the function, and you just didn't update the doc with the code.

Other than that, providing a docstring is good. I'm not convinced by the idea of having sample values in the docstring instead of (as well as would be fine) variable names and explanations. (Of course it did happen to be very helpful for my profiling to have an input I knew should work!)

---

Python style is pretty well defined, although deviations are allowed if they're applied consistently. Function names should be lowercase and snake_case, as indeed you have with everything except CycloPepSeq. Likewise variable names, which you normally do and occasionally fall back to PrefixMass and such.

It's worth getting hold of Pylint or similar to help enforce conventions like this, just to make reading the code more natural for other python programmers.

That said, don't be a slave to it. For example it also objects to the long line defining AAmass. Given what that line is doing, it's plausibly better off stretching off the screen. After all it isn't really code in the sense that people have to comprehend it, beyond the vague gist of "lots of data".

---

One potentially more significant issue pylint discovered is the redefinition of score as a variable within score as a function. It hasn't done any harm in this instance because scoping rules protected it, but it can cause quite confusing bugs.

---

Finally, tests would be good.

Certainly once you have something that works, you shouldn't dream of doing an extensive refactor to optimise it until such time as you have unit and integration tests to keep you on the right track. Ideally they'd cover a range of possible inputs, including short ones (lengths 0, 1, 2) and a spread of long ones, as well as edge cases like input lists with lots of duplicate vaues.

---

Testing against the large input now provided, a single run of the original program took 40.4 seconds on my machine as written, of which score represented 22.8 seconds and remove within score a further 10.1 seconds. Applying the changes above, but replacing set with collections.Counter and return len(idealspectrum & spectrum) with return sum((spectrum_counter & idealspectrum).values()) gives the same output for that large example and completes in 11.7 seconds.