Mass spectrometry

5 problems from Rosalind — Bioinformatics Stronghold. Press Run on any block to execute it in your browser.

# Rosalind: SPEC — Inferring Protein from Spectrum
# https://rosalind.info/problems/spec/
#
# Given: A list L of n masses forming the prefix spectrum of a protein.
# Return: A protein string whose prefix spectrum is L.

let spectrum = [3524.8542, 3710.9335, 3841.974, 3970.0326, 4041.0697]

let monoisotopic = {
    A: 71.03711,  C: 103.00919, D: 115.02694, E: 129.04259, F: 147.06841,
    G: 57.02146,  H: 137.05891, I: 113.08406, K: 128.09496, L: 113.08406,
    M: 131.04049, N: 114.04293, P: 97.05276,  Q: 128.05858, R: 156.10111,
    S: 87.03203,  T: 101.04768, V: 99.06841,  W: 186.07931, Y: 163.06333
}

# Consecutive masses in a prefix spectrum differ by exactly one residue, so
# each gap identifies the residue nearest it. Leucine and isoleucine share a
# mass; the ordering of `keys` decides which name is reported.
fn residue_for(gap) {
    let names = keys(monoisotopic)
    let scored = names |> map(|name| { name: name, error: abs(monoisotopic[name] - gap) })
    let best = scored |> sort_by(|e| e.error) |> first()
    best.name
}

let peptide = range(1, len(spectrum))
    |> map(|i| residue_for(spectrum[i] - spectrum[i - 1]))
    |> join("")

println("Result:   " + peptide)
println("Gaps:     " + (range(1, len(spectrum)) |> map(|i| str(round(spectrum[i] - spectrum[i - 1], 4))) |> join(" ")))
println("Expected: WMQA — 186.0793 W, 131.0405 M, 128.0586 Q, 71.0371 A")

fn test_spec_protein_from_spectrum() {
    assert peptide == "WMQA", "SPEC: got " + peptide
    # Every reported residue must actually reproduce its gap, which is the
    # property the answer rests on.
    let i = 1
    while i < len(spectrum) {
        let gap = spectrum[i] - spectrum[i - 1]
        let residue = substr(peptide, i - 1, 1)
        assert abs(monoisotopic[residue] - gap) < 0.001, "SPEC: residue " + residue + " does not match its gap"
        i = i + 1
    }
}

# Rosalind: CONV — Comparing Spectra with the Spectral Convolution
# https://rosalind.info/problems/conv/
#
# Given: Two multisets of masses S1 and S2.
# Return: The largest multiplicity of the spectral convolution S1 (-) S2,
# followed by the value achieving it.

let s1 = [186.07931, 287.12699, 548.20532, 580.18077, 681.22845, 706.27446, 782.27613, 968.35544, 968.35544]
let s2 = [101.04768, 158.06914, 202.09536, 318.09979, 419.14747, 463.17369, 507.19992, 536.21545,
          597.25729, 618.28871, 664.27596, 682.25123, 785.29975, 787.28104, 803.29542, 819.28958,
          819.28958, 891.35053, 924.36613, 1069.44506]

# The convolution is every pairwise difference. Masses are compared at five
# decimal places: they come from instrument readings, so exact float equality
# would split values that are meant to be the same.
let differences = s1 |> flat_map(|a| s2 |> map(|b| round(a - b, 5)))

let tallied = differences |> unique() |> map(|d| {
    { value: d, count: differences |> count_if(|x| x == d) }
})

let best = tallied |> sort_by(|e| e.count) |> reverse() |> first()

let count_85 = differences |> count_if(|d| d == 85.03163)

println("Result:   " + str(best.count))
println("          " + str(best.value))
println("")
println("Both spectra repeat a value — 968.35544 twice in S1, 819.28958 twice")
println("in S2 — and their difference is 149.06586, so that pairing occurs")
println("2 x 2 = 4 times. Multiplicity is over the multiset, so it wins.")
println("")
println("85.03163 occurs " + str(count_85) + " times:")
println("  186.07931-101.04768, 287.12699-202.09536, 548.20532-463.17369")

fn test_conv_spectral_convolution() {
    # The duplicated masses pair four ways.
    assert best.count == 4, "CONV: multiplicity " + str(best.count)
    assert best.value == 149.06586, "CONV: value " + str(best.value)
    # And the three hand-checked differences are all present.
    assert count_85 == 3, "CONV: 85.03163 occurred " + str(count_85) + " times"
    # Every difference must be reproducible from some pair of inputs.
    assert len(differences) == len(s1) * len(s2), "CONV: convolution size is wrong"
}

# Rosalind: PRSM — Matching a Spectrum to a Protein
# https://rosalind.info/problems/prsm/
#
# Given: A collection of protein strings and a multiset R of masses.
# Return: The largest multiplicity of any protein's complete spectrum convolved
# with R, and a protein achieving it.

let proteins = ["GSDMQS", "VWICN", "IASMQS", "PVSMGAD"]
let spectrum = [445.17838, 115.02694, 186.07931, 314.13789, 317.1198, 215.09061]

let monoisotopic = {
    A: 71.03711,  C: 103.00919, D: 115.02694, E: 129.04259, F: 147.06841,
    G: 57.02146,  H: 137.05891, I: 113.08406, K: 128.09496, L: 113.08406,
    M: 131.04049, N: 114.04293, P: 97.05276,  Q: 128.05858, R: 156.10111,
    S: 87.03203,  T: 101.04768, V: 99.06841,  W: 186.07931, Y: 163.06333
}

fn mass_of(peptide) {
    range(0, len(peptide)) |> reduce(|acc, i| acc + monoisotopic[substr(peptide, i, 1)], 0.0)
}

# The complete spectrum is every prefix and every suffix mass — the fragments a
# spectrometer would see if the peptide broke at each bond in turn.
fn complete_spectrum(peptide) {
    let prefixes = range(1, len(peptide)) |> map(|k| mass_of(substr(peptide, 0, k)))
    let suffixes = range(1, len(peptide)) |> map(|k| mass_of(substr(peptide, len(peptide) - k, k)))
    concat(prefixes, suffixes)
}

# Convolving the two multisets, the winning shift is the parent mass offset;
# its multiplicity is how many fragments line up.
fn best_multiplicity(peptide) {
    let own = complete_spectrum(peptide)
    let shifts = spectrum |> flat_map(|r| own |> map(|m| round(r - m, 5)))
    let tallies = shifts |> unique() |> map(|d| shifts |> count_if(|x| x == d))
    max(tallies)
}

let scored = proteins |> map(|p| { protein: p, score: best_multiplicity(p) })
let best = scored |> sort_by(|e| e.score) |> reverse() |> first()

println("Scores:")
scored |> each(|e| println("  " + e.protein + ": " + str(e.score)))
println("Result:   " + str(best.score))
println("          " + best.protein)
println("Expected: 3 and IASMQS")
println("")
println("Note: GSDMQS also reaches 3, so the maximum is shared. Rosalind accepts")
println("any protein achieving it; which one is reported here depends on the")
println("order of the input rather than on anything meaningful.")

let maximal = scored |> filter(|e| e.score == best.score) |> map(|e| e.protein)

fn test_prsm_spectrum_to_protein() {
    assert best.score == 3, "PRSM: multiplicity " + str(best.score)
    # The tie is real, so require the winner to be one of the maximal proteins
    # rather than one particular string.
    assert maximal |> contains(best.protein), "PRSM: winner is not maximal"
    assert maximal |> contains("IASMQS"), "PRSM: IASMQS is not among the maximal proteins"
    # A complete spectrum holds two fragments per bond.
    assert len(complete_spectrum("IASMQS")) == 2 * (len("IASMQS") - 1), "PRSM: spectrum size is wrong"
}

# Rosalind: SGRA — Using the Spectrum Graph to Infer Peptides
# https://rosalind.info/problems/sgra/
#
# Given: A list of positive real numbers (a spectrum).
# Return: The longest protein string matching the spectrum graph.

let spectrum = [
    3524.8542, 3623.5245, 3710.9335, 3841.974, 3929.00603,
    3970.0326, 4026.05879, 4057.0646, 4083.08025,
]

# Monoisotopic masses, not the integer table the cyclopeptide problems use.
# Real instrument readings carry decimals, so a difference is matched within a
# tolerance rather than exactly — 1e-4 here, which is far tighter than the gap
# between any two residue masses and far looser than float noise.
let monoisotopic = {
    A: 71.03711,  C: 103.00919, D: 115.02694, E: 129.04259, F: 147.06841,
    G: 57.02146,  H: 137.05891, I: 113.08406, K: 128.09496, L: 113.08406,
    M: 131.04049, N: 114.04293, P: 97.05276,  Q: 128.05858, R: 156.10111,
    S: 87.03203,  T: 101.04768, V: 99.06841,  W: 186.07931, Y: 163.06333
}
let tolerance = 0.0001

# Where L and I collide exactly, and K and Q nearly do, take the first — a
# spectrum cannot separate them, so naming either is equally correct.
let residues = ["G", "A", "S", "P", "V", "T", "C", "L", "N", "D",
                "Q", "K", "E", "M", "H", "F", "R", "Y", "W"]

fn residue_between(lighter, heavier, table, letters, slack) {
    let gap = heavier - lighter
    let hits = letters |> filter(|letter| abs(table[letter] - gap) < slack)
    if len(hits) > 0 then hits[0] else ""
}

# The graph is a DAG — masses only increase — so the longest path is found by
# scanning in sorted order and never revisiting a node. Searching the paths
# themselves would be exponential; here each node is settled once, from the best
# way of reaching it.
let ordered = sort(spectrum)
let best = ordered |> map(|_| "")
# Where each best path began, so the answer's mass can be checked against the
# span it was actually read across. The path skips nodes, so this cannot be
# recovered by counting backwards from the end.
let origin = range(0, len(ordered))

for j in range(0, len(ordered)) {
    for i in range(0, j) {
        let letter = residue_between(ordered[i], ordered[j], monoisotopic, residues, tolerance)
        if letter != "" and len(best[i]) + 1 > len(best[j]) {
            best[j] = best[i] + letter
            origin[j] = origin[i]
        }
    }
}

let finish = argmax(best |> map(|s| len(s)))
let answer = best[finish]

println("Result:   " + answer)
println("Expected: WMSPG")

fn test_sgra_spectrum_graph() {
    assert answer == "WMSPG", "SGRA: got " + answer
    # Every residue of the answer must be a real gap between two spectrum values.
    assert len(answer) == 5, "SGRA: expected a 5-residue peptide"
    # The peptide's own mass is the span it was read across.
    let peptide_mass = chars(answer) |> map(|c| monoisotopic[c]) |> sum()
    let span = ordered[finish] - ordered[origin[finish]]
    assert abs(peptide_mass - span) < len(answer) * tolerance,
        "SGRA: the peptide weighs " + str(round(peptide_mass, 4))
            + " but spans " + str(round(span, 4))
    # Nothing longer exists in the graph.
    for candidate in best {
        assert len(candidate) <= len(answer), "SGRA: a longer path was missed"
    }
}

# Rosalind: FULL — Inferring Peptide from Full Spectrum
# https://rosalind.info/problems/full/
#
# Given: 2n+3 positive reals — a parent mass, then the b-ions and y-ions of a
# peptide of length n, in no particular order.
# Return: A protein string of length n consistent with them.

let parent_mass = 1988.21104821
let ions = [
    610.391039105, 738.485999105, 766.492149105, 863.544909105,
    867.528589105, 992.587499105, 995.623549105, 1120.6824591,
    1124.6661391, 1221.7188991, 1249.7250491, 1377.8200091,
]

let monoisotopic = {
    A: 71.03711,  C: 103.00919, D: 115.02694, E: 129.04259, F: 147.06841,
    G: 57.02146,  H: 137.05891, I: 113.08406, K: 128.09496, L: 113.08406,
    M: 131.04049, N: 114.04293, P: 97.05276,  Q: 128.05858, R: 156.10111,
    S: 87.03203,  T: 101.04768, V: 99.06841,  W: 186.07931, Y: 163.06333
}
let residues = ["G", "A", "S", "P", "V", "T", "C", "L", "N", "D",
                "Q", "K", "E", "M", "H", "F", "R", "Y", "W"]
let tolerance = 0.001

# A peptide fragments from both ends at once. A b-ion is a prefix, a y-ion the
# matching suffix, so the two always sum to the parent mass — and the list mixes
# them with no labels saying which is which.
#
# The trick is that it does not matter. Walk upwards from the lightest ion; the
# next ion differing from it by a residue mass is the next prefix, whichever kind
# it nominally is. Taking an ion also spends its complement, since a prefix and
# its suffix are one fragmentation event and cannot both extend the chain.
let ordered = sort(ions)
let used = ordered |> map(|_| false)

fn residue_between(lighter, heavier, table, letters, slack) {
    let gap = heavier - lighter
    let hits = letters |> filter(|letter| abs(table[letter] - gap) < slack)
    if len(hits) > 0 then hits[0] else ""
}

fn index_of_mass(masses, wanted, slack) {
    let hits = range(0, len(masses)) |> filter(|i| abs(masses[i] - wanted) < slack)
    if len(hits) > 0 then hits[0] else 0 - 1
}

fn spend(marks, masses, at, total, slack) {
    marks[at] = true
    let partner = index_of_mass(masses, total - masses[at], slack)
    if partner >= 0 { marks[partner] = true }
    marks
}

let peptide = ""
let current = 0
used = spend(used, ordered, current, parent_mass, tolerance)

let searching = true
while searching {
    let onward = range(current + 1, len(ordered))
        |> filter(|j| used[j] == false
                  and residue_between(ordered[current], ordered[j],
                                      monoisotopic, residues, tolerance) != "")
    if len(onward) == 0 {
        searching = false
    } else {
        let next_index = onward[0]
        peptide = peptide + residue_between(ordered[current], ordered[next_index],
                                            monoisotopic, residues, tolerance)
        used = spend(used, ordered, next_index, parent_mass, tolerance)
        current = next_index
    }
}

println("Result:   " + peptide)
println("Expected: KEKEP")

fn test_full_inferring_peptide_from_full_spectrum() {
    assert peptide == "KEKEP", "FULL: got " + peptide
    # n is fixed by the input's size, and the answer has to match it.
    let n = (len(ions) + 1 - 3) / 2
    assert len(peptide) == n, "FULL: expected " + str(n) + " residues"
    # Every ion pairs with another summing to the parent mass — the property that
    # makes b-ions and y-ions indistinguishable here, and the reason taking one
    # spends the other.
    for mass in ordered {
        let partner = index_of_mass(ordered, parent_mass - mass, tolerance)
        assert partner >= 0, "FULL: " + str(round(mass, 4)) + " has no complement"
    }
    assert len(ordered) % 2 == 0, "FULL: the ions must pair up"
}