nx_peptide_denovo.nx source
↩ module page · 162 lines · 6504 B
1// nx_peptide_denovo.nx -- SOTA: DE NOVO peptide sequencing. Reconstruct the
2// amino-acid sequence from a fragment-ion spectrum ALONE, with NO candidate
3// and NO database. This is the hardest problem in proteomics: read the
4// residues directly off the mass GAPS in the b-ion ladder.
5//
6// THE PRINCIPLE, EXACT. Consecutive b-ions differ by exactly one residue's
7// mass: b_i - b_(i-1) = residue i. The first residue is b_1 - proton; the
8// last is recovered from the precursor mass. So a sorted b-ion ladder plus
9// the intact mass spells the sequence, gap by gap -- match each gap to the
10// nearest residue mass and read the letter. nx_peptide_identify needed a
11// candidate to score against; this needs nothing but the spectrum.
12//
13// THE HONEST AMBIGUITIES, ENCODED NOT HIDDEN. Two residue pairs cannot be
14// told apart by mass and de novo sequencing is famous for exactly this:
15// Leu / Ile -- IDENTICAL mass. De novo cannot distinguish them; this
16// returns a canonical 'L' and flags the position.
17// Lys / Gln -- differ by 0.036 Da; separable only at high mass accuracy.
18// A de novo tool that pretended to resolve Leu/Ile would be lying about
19// physics. This one reports the ambiguity rather than guessing.
20//
21// ANCHORED. Residue masses are the anchored table two independent methods
22// agree on; the gap arithmetic is exact. The gate liar-kills it by ROUND
23// TRIP: take a known peptide, generate its real b-ion ladder, de novo
24// sequence it back, and recover the original sequence (modulo Leu/Ile). If a
25// residue-mass gap were wrong, a letter would come back wrong.
26//
27// All INTEGER, _q4 = x10^4 Da.
28//
29// Grounding (cited; researcher-groundable):
30// dancik_1999_de_novo_peptide_sequencing
31// biemann_1990_fragment_ion_nomenclature
32// leu_ile_isobaric_lys_gln_near_isobaric
33//
34// genealogy_id: peptide_chemistry + analytical_sota + liar_killer
35
36import "nx_syscalls.nx"
37import "nx_peptide.nx"
38
39const DENOVO_UNKNOWN: i64 = 0 - 1 // no residue matches a gap
40const DENOVO_MAX_LEN: i64 = 128
41
42// ===== Residue from a mass gap ========================================
43//
44// Which residue has monoisotopic mass within tol of delta? Table order means
45// Leu (id 9) is checked before Ile (id 7)? No -- we WANT a canonical choice
46// for the isobaric pair, so check the canonical FIRST. Returns the aa id, or
47// DENOVO_UNKNOWN. Leu/Ile share a mass; this returns PEP_L canonically.
48
49func denovo_residue(delta_q4: i64, tol_q4: i64) -> i64 {
50 var aa: i64 = 0
51 var best: i64 = DENOVO_UNKNOWN
52 var best_d: i64 = tol_q4 + 1
53 // Check Leu first so the isobaric Leu/Ile pair resolves to a canonical L.
54 let lm: i64 = pep_residue_mono_q4(PEP_L)
55 var dl: i64 = delta_q4 - lm
56 if dl < 0 { dl = 0 - dl }
57 if dl <= tol_q4 { best = PEP_L; best_d = dl }
58 while aa < PEP_N_AA {
59 if aa != PEP_I {
60 if aa != PEP_L {
61 let m: i64 = pep_residue_mono_q4(aa)
62 var d: i64 = delta_q4 - m
63 if d < 0 { d = 0 - d }
64 if d <= tol_q4 { if d < best_d { best = aa; best_d = d } }
65 }
66 }
67 aa = aa + 1
68 }
69 return best
70}
71
72// Is a residue at this position ambiguous (Leu/Ile isobaric)?
73func denovo_is_leu_ile(delta_q4: i64, tol_q4: i64) -> i64 {
74 let lm: i64 = pep_residue_mono_q4(PEP_L)
75 var d: i64 = delta_q4 - lm
76 if d < 0 { d = 0 - d }
77 if d <= tol_q4 { return 1 }
78 return 0
79}
80
81// aa id -> 1-letter ASCII code (inverse of pep_aa_from_char).
82func denovo_char(aa: i64) -> i64 {
83 if aa == PEP_A { return 65 }
84 if aa == PEP_C { return 67 }
85 if aa == PEP_D { return 68 }
86 if aa == PEP_E { return 69 }
87 if aa == PEP_F { return 70 }
88 if aa == PEP_G { return 71 }
89 if aa == PEP_H { return 72 }
90 if aa == PEP_I { return 73 }
91 if aa == PEP_K { return 75 }
92 if aa == PEP_L { return 76 }
93 if aa == PEP_M { return 77 }
94 if aa == PEP_N { return 78 }
95 if aa == PEP_P { return 80 }
96 if aa == PEP_Q { return 81 }
97 if aa == PEP_R { return 82 }
98 if aa == PEP_S { return 83 }
99 if aa == PEP_T { return 84 }
100 if aa == PEP_V { return 86 }
101 if aa == PEP_W { return 87 }
102 if aa == PEP_Y { return 89 }
103 return 88 // 'X' -- unknown
104}
105
106// ===== De novo sequence ===============================================
107//
108// From a SORTED b-ion ladder (b_1 < b_2 < ... < b_(n-1)) and the precursor
109// neutral mass, reconstruct the n-residue sequence into out (NUL-terminated).
110// Returns the length, or DENOVO_UNKNOWN if any gap has no residue match.
111//
112// residue 1 = b_1 - proton
113// residue i = b_i - b_(i-1) (2 <= i <= n-1)
114// residue n = precursor - water - (b_(n-1) - proton)
115
116func denovo_sequence(b_ions: *i64, n_b: i64, precursor_q4: i64, tol_q4: i64, out: *u8) -> i64 {
117 if n_b < 1 { return DENOVO_UNKNOWN }
118 if n_b + 1 > DENOVO_MAX_LEN { return DENOVO_UNKNOWN }
119 var pos: i64 = 0
120 // residue 1 from b_1
121 let r1: i64 = b_ions[0] - PEP_PROTON_Q4
122 let aa1: i64 = denovo_residue(r1, tol_q4)
123 if aa1 == DENOVO_UNKNOWN { return DENOVO_UNKNOWN }
124 out[pos] = denovo_char(aa1) as u8
125 pos = pos + 1
126 // residues 2..n-1 from consecutive gaps
127 var i: i64 = 1
128 while i < n_b {
129 let delta: i64 = b_ions[i] - b_ions[i - 1]
130 let aa: i64 = denovo_residue(delta, tol_q4)
131 if aa == DENOVO_UNKNOWN { return DENOVO_UNKNOWN }
132 out[pos] = denovo_char(aa) as u8
133 pos = pos + 1
134 i = i + 1
135 }
136 // last residue from the precursor mass
137 let last_gap: i64 = precursor_q4 - PEP_WATER_MONO_Q4 - (b_ions[n_b - 1] - PEP_PROTON_Q4)
138 let aal: i64 = denovo_residue(last_gap, tol_q4)
139 if aal == DENOVO_UNKNOWN { return DENOVO_UNKNOWN }
140 out[pos] = denovo_char(aal) as u8
141 pos = pos + 1
142 out[pos] = 0 as u8
143 return pos
144}
145
146// How many positions are Leu/Ile-ambiguous in a b-ion ladder (the honest
147// uncertainty count a de novo call must report).
148func denovo_ambiguous_count(b_ions: *i64, n_b: i64, precursor_q4: i64, tol_q4: i64) -> i64 {
149 if n_b < 1 { return 0 }
150 var amb: i64 = 0
151 let r1: i64 = b_ions[0] - PEP_PROTON_Q4
152 amb = amb + denovo_is_leu_ile(r1, tol_q4)
153 var i: i64 = 1
154 while i < n_b {
155 let delta: i64 = b_ions[i] - b_ions[i - 1]
156 amb = amb + denovo_is_leu_ile(delta, tol_q4)
157 i = i + 1
158 }
159 let last_gap: i64 = precursor_q4 - PEP_WATER_MONO_Q4 - (b_ions[n_b - 1] - PEP_PROTON_Q4)
160 amb = amb + denovo_is_leu_ile(last_gap, tol_q4)
161 return amb
162}