nx_peptide_msms.nx source
↩ module page · 220 lines · 8083 B
1// nx_peptide_msms.nx -- SOTA ANALYTICAL CAPSTONE: peptide identification by
2// MS/MS spectral matching. Given an observed fragment-ion peak list, it
3// generates the theoretical b/y ladder for a candidate sequence, matches
4// observed against theoretical within a tolerance, and returns an identity
5// score, the sequence coverage, and -- across a set of candidates -- which
6// peptide the spectrum IS. This is the literal core of how a mass
7// spectrometer identifies a peptide.
8//
9// WHY IT IS STATE OF THE ART. A single-mass match ("a peak near 1046") is a
10// guess; a fragment-spectrum match is an identification. Backbone cleavage
11// gives the b (N-terminal) and y (C-terminal) ion ladders, and a real ID
12// scores how much of that predicted ladder the instrument actually saw, and
13// how much of the backbone that covers. That is what nx_peptide_msms computes.
14//
15// ANCHORED, NOT ASSERTED. The b/y ion masses come from the residue table that
16// two independent methods (residue-sum and elemental formula) agree on and
17// that matches published values; the scoring is exact integer arithmetic. The
18// gate liar-kills it: a spectrum matched against its OWN sequence scores a
19// perfect 1000, and against a DIFFERENT peptide scores low -- if the matcher
20// could not tell them apart it would be useless, and the negative control
21// proves it can.
22//
23// The b_i + y_(n-i) = M + 2*proton invariant (proven in nx_peptide) means the
24// two ladders are not independent evidence of the same thing twice; coverage
25// counts a backbone SITE as seen if EITHER its b or its complementary y ion
26// appears, which is the honest way to score.
27//
28// All INTEGER, _q4 = x10^4 Da matching nx_peptide.
29//
30// Grounding (cited; researcher-groundable):
31// biemann_1990_peptide_fragment_ion_nomenclature
32// eng_1994_sequest_cross_correlation_spectral_match
33// published_peptide_ms_ms_fragment_ladders
34//
35// genealogy_id: peptide_chemistry + analytical_sota + liar_killer
36
37import "nx_syscalls.nx"
38import "nx_peptide.nx"
39
40const MSMS_MAX_IONS: i64 = 512 // 2*(n-1); a 256-mer is well past any real peptide
41
42// ===== Theoretical spectrum ===========================================
43//
44// Emit the full b/y ladder (b1..b(n-1), y1..y(n-1)) into out_mz. Returns the
45// ion count, or MSMS refusal (-1) on a bad sequence / overflow.
46
47const MSMS_REFUSED: i64 = 0 - 1
48
49func msms_theoretical_spectrum(seq: *u8, out_mz: *i64, cap: i64) -> i64 {
50 if pep_seq_valid(seq) != 1 { return MSMS_REFUSED }
51 let n: i64 = pep_len(seq)
52 if n < 2 { return MSMS_REFUSED }
53 var count: i64 = 0
54 var i: i64 = 1
55 while i < n {
56 let b: i64 = pep_b_ion_q4(seq, i)
57 let y: i64 = pep_y_ion_q4(seq, i)
58 if b != PEP_INVALID {
59 if count < cap { out_mz[count] = b; count = count + 1 }
60 }
61 if y != PEP_INVALID {
62 if count < cap { out_mz[count] = y; count = count + 1 }
63 }
64 i = i + 1
65 }
66 return count
67}
68
69// Number of DISTINCT theoretical ions for a length-n peptide = 2*(n-1).
70func msms_theoretical_count(seq: *u8) -> i64 {
71 let n: i64 = pep_len(seq)
72 if n < 2 { return 0 }
73 return 2 * (n - 1)
74}
75
76// ===== Matching =======================================================
77
78// Is target within tol of ANY peak in the observed list?
79func msms_peak_present(target: i64, obs: *i64, n_obs: i64, tol_q4: i64) -> i64 {
80 var k: i64 = 0
81 var hit: i64 = 0
82 while k < n_obs {
83 if hit == 0 {
84 var d: i64 = obs[k] - target
85 if d < 0 { d = 0 - d }
86 if d <= tol_q4 { hit = 1 }
87 }
88 k = k + 1
89 }
90 return hit
91}
92
93// How many theoretical b/y ions of `seq` appear in the observed list.
94func msms_matched_ions(seq: *u8, obs: *i64, n_obs: i64, tol_q4: i64) -> i64 {
95 if pep_seq_valid(seq) != 1 { return 0 }
96 let n: i64 = pep_len(seq)
97 if n < 2 { return 0 }
98 var matched: i64 = 0
99 var i: i64 = 1
100 while i < n {
101 let b: i64 = pep_b_ion_q4(seq, i)
102 let y: i64 = pep_y_ion_q4(seq, i)
103 if b != PEP_INVALID { matched = matched + msms_peak_present(b, obs, n_obs, tol_q4) }
104 if y != PEP_INVALID { matched = matched + msms_peak_present(y, obs, n_obs, tol_q4) }
105 i = i + 1
106 }
107 return matched
108}
109
110// Identity score: matched theoretical ions / total theoretical ions, permil.
111// A perfect self-match is 1000; an unrelated spectrum is near 0.
112func msms_score_permil(seq: *u8, obs: *i64, n_obs: i64, tol_q4: i64) -> i64 {
113 let total: i64 = msms_theoretical_count(seq)
114 if total <= 0 { return 0 }
115 let matched: i64 = msms_matched_ions(seq, obs, n_obs, tol_q4)
116 return matched * 1000 / total
117}
118
119// ===== Sequence coverage ==============================================
120//
121// A backbone cleavage SITE i is covered if its b_i OR its complementary
122// y_(n-i) ion was observed. Coverage = covered sites / (n-1) backbone bonds,
123// permil. Counting the site once (not b and y separately) is honest, since
124// b_i and y_(n-i) are complementary evidence of the SAME cleavage.
125
126func msms_coverage_permil(seq: *u8, obs: *i64, n_obs: i64, tol_q4: i64) -> i64 {
127 if pep_seq_valid(seq) != 1 { return 0 }
128 let n: i64 = pep_len(seq)
129 if n < 2 { return 0 }
130 let sites: i64 = n - 1
131 var covered: i64 = 0
132 var i: i64 = 1
133 while i < n {
134 let b: i64 = pep_b_ion_q4(seq, i)
135 let y: i64 = pep_y_ion_q4(seq, n - i)
136 var seen: i64 = 0
137 if b != PEP_INVALID {
138 if msms_peak_present(b, obs, n_obs, tol_q4) == 1 { seen = 1 }
139 }
140 if y != PEP_INVALID {
141 if msms_peak_present(y, obs, n_obs, tol_q4) == 1 { seen = 1 }
142 }
143 covered = covered + seen
144 i = i + 1
145 }
146 return covered * 1000 / sites
147}
148
149// ===== Identification across candidates ================================
150//
151// A minimal library search: score a spectrum against K candidate sequences
152// and return the index of the best. Ties go to the earlier candidate. This is
153// how a spectrum becomes an identity.
154
155struct NxCandidates {
156 seqs: *i64, // array of *u8 (sequence pointers)
157 n: i64,
158}
159
160func nx_candidates_new(cap: i64) -> *NxCandidates {
161 let c: *NxCandidates = (sys_mmap(24)) as *NxCandidates
162 c.seqs = sys_mmap(cap * 8)
163 c.n = 0
164 return c
165}
166
167func msms_add_candidate(c: *NxCandidates, seq: *u8) -> i64 {
168 let s: *i64 = c.seqs
169 s[c.n] = seq as i64
170 c.n = c.n + 1
171 return c.n
172}
173
174func msms_best_candidate(c: *NxCandidates, obs: *i64, n_obs: i64, tol_q4: i64) -> i64 {
175 let s: *i64 = c.seqs
176 var best: i64 = 0 - 1
177 var bestscore: i64 = 0 - 1
178 var i: i64 = 0
179 while i < c.n {
180 let seq: *u8 = s[i] as *u8
181 let sc: i64 = msms_score_permil(seq, obs, n_obs, tol_q4)
182 if sc > bestscore { bestscore = sc; best = i }
183 i = i + 1
184 }
185 return best
186}
187
188func msms_best_score_permil(c: *NxCandidates, obs: *i64, n_obs: i64, tol_q4: i64) -> i64 {
189 let idx: i64 = msms_best_candidate(c, obs, n_obs, tol_q4)
190 if idx < 0 { return 0 }
191 let s: *i64 = c.seqs
192 let seq: *u8 = s[idx] as *u8
193 return msms_score_permil(seq, obs, n_obs, tol_q4)
194}
195
196// Identity confidence: an ID is trustworthy only if the best score clears a
197// floor AND beats the runner-up by a margin (a spectrum that matches two
198// candidates equally has identified nothing). Returns 1 iff confident.
199const MSMS_SCORE_FLOOR: i64 = 500 // >=50% of the ladder seen
200const MSMS_MARGIN_PERMIL: i64 = 200 // best beats 2nd by >=20 points
201
202func msms_confident_id(c: *NxCandidates, obs: *i64, n_obs: i64, tol_q4: i64) -> i64 {
203 if c.n < 1 { return 0 }
204 let s: *i64 = c.seqs
205 var best: i64 = 0
206 var second: i64 = 0
207 var i: i64 = 0
208 while i < c.n {
209 let seq: *u8 = s[i] as *u8
210 let sc: i64 = msms_score_permil(seq, obs, n_obs, tol_q4)
211 if sc > best { second = best; best = sc }
212 if sc <= best {
213 if sc > second { if sc < best { second = sc } }
214 }
215 i = i + 1
216 }
217 if best < MSMS_SCORE_FLOOR { return 0 }
218 if best - second < MSMS_MARGIN_PERMIL { return 0 }
219 return 1
220}