nx_peptide_msms_test.nx source
↩ module page · 145 lines · 8941 B
1// nx_peptide_msms_test.nx -- gate for the MS/MS identification capstone. The
2// liar-kill: a spectrum scored against its OWN sequence is a perfect 1000, and
3// against a DIFFERENT peptide is low (T1/T2). If the matcher could not
4// discriminate it would be worthless; the negative control proves it can.
5// expect_exit: 0 license_tier: ORIGINAL
6import "nx_syscalls.nx"
7import "nx_peptide.nx"
8import "nx_peptide_msms.nx"
9
10func t_puts(s: *u8) -> i64 { var n: i64 = 0; while s[n] != (0 as u8) { n = n + 1 } sys_write(1, s, n); return 0 }
11func t_putn(v: i64) -> i64 { let bb: *u8 = sys_mmap(28); var m: i64 = v; if m < 0 { m = 0 - m; sys_write(1, "-" as *u8, 1) } let t: *u8 = sys_mmap(28); var k: i64 = 0; if m == 0 { t[0] = 48 as u8; k = 1 } while m > 0 { t[k] = (48 + (m % 10)) as u8; m = m / 10; k = k + 1 } var i: i64 = 0; while i < k { bb[i] = t[k - 1 - i]; i = i + 1 } sys_write(1, bb, k); return 0 }
12
13func main() -> i64 {
14 var pass: i64 = 0
15 var total: i64 = 0
16
17 let ang: *u8 = "DRVYIHPF" as *u8
18 let bk: *u8 = "RPPGFSPFR" as *u8
19 let tol: i64 = 100 // 0.01 Da
20
21 // Generate angiotensin's own theoretical spectrum -> use it as "observed".
22 let ang_spec: *i64 = sys_mmap(MSMS_MAX_IONS * 8)
23 let ang_n: i64 = msms_theoretical_spectrum(ang, ang_spec, MSMS_MAX_IONS)
24
25 // --- T1 PERFECT SELF-MATCH: angiotensin vs its own spectrum = 1000 ---
26 total = total + 1
27 let self_score: i64 = msms_score_permil(ang, ang_spec, ang_n, tol)
28 let self_cov: i64 = msms_coverage_permil(ang, ang_spec, ang_n, tol)
29 t_puts("T1 angiotensin vs OWN spectrum: score=" as *u8); t_putn(self_score); t_puts(" coverage=" as *u8); t_putn(self_cov); t_puts(" want 1000/1000: " as *u8)
30 var ok1: i64 = 1
31 if self_score != 1000 { ok1 = 0 }
32 if self_cov != 1000 { ok1 = 0 }
33 if ok1 == 1 { pass = pass + 1; t_puts("PASS\n" as *u8) } else { t_puts("FAIL\n" as *u8) }
34
35 // --- T2 DISCRIMINATION (neg-control): bradykinin scored against
36 // angiotensin's spectrum must be LOW ---
37 total = total + 1
38 let cross: i64 = msms_score_permil(bk, ang_spec, ang_n, tol)
39 t_puts("T2 bradykinin vs angiotensin spectrum: score=" as *u8); t_putn(cross); t_puts(" permil (must be low, <300): " as *u8)
40 if cross < 300 { pass = pass + 1; t_puts("PASS\n" as *u8) } else { t_puts("FAIL\n" as *u8) }
41
42 // --- T3 the theoretical spectrum has 2*(n-1) ions for an 8-mer = 14 ---
43 total = total + 1
44 let tc: i64 = msms_theoretical_count(ang)
45 t_puts("T3 angiotensin (8-mer) theoretical ion count=" as *u8); t_putn(ang_n); t_puts(" formula 2*(n-1)=" as *u8); t_putn(tc); t_puts(" want 14/14: " as *u8)
46 var ok3: i64 = 1
47 if ang_n != 14 { ok3 = 0 }
48 if tc != 14 { ok3 = 0 }
49 if ok3 == 1 { pass = pass + 1; t_puts("PASS\n" as *u8) } else { t_puts("FAIL\n" as *u8) }
50
51 // --- T4 a PARTIAL spectrum (only the first half of the ions) scores
52 // lower but still non-trivially -- a real spectrum is never complete ---
53 total = total + 1
54 let half_n: i64 = ang_n / 2
55 let half_score: i64 = msms_score_permil(ang, ang_spec, half_n, tol)
56 t_puts("T4 angiotensin vs HALF its spectrum: score=" as *u8); t_putn(half_score); t_puts(" permil (partial, 300..700): " as *u8)
57 var ok4: i64 = 1
58 if half_score < 300 { ok4 = 0 }
59 if half_score > 700 { ok4 = 0 }
60 if half_score >= self_score { ok4 = 0 }
61 if ok4 == 1 { pass = pass + 1; t_puts("PASS\n" as *u8) } else { t_puts("FAIL\n" as *u8) }
62
63 // --- T5 IDENTIFICATION: a library search over candidates picks the
64 // right peptide from an observed spectrum ---
65 total = total + 1
66 let lib: *NxCandidates = nx_candidates_new(8)
67 msms_add_candidate(lib, bk) // 0 bradykinin
68 msms_add_candidate(lib, "HWAWFK" as *u8) // 1 GHRP-6 backbone
69 msms_add_candidate(lib, ang) // 2 angiotensin (the truth)
70 msms_add_candidate(lib, "LKKTETQ" as *u8)// 3 TB-500 fragment
71 let best: i64 = msms_best_candidate(lib, ang_spec, ang_n, tol)
72 let best_score: i64 = msms_best_score_permil(lib, ang_spec, ang_n, tol)
73 t_puts("T5 library search on angiotensin spectrum -> candidate " as *u8); t_putn(best); t_puts(" (want 2=angiotensin) score=" as *u8); t_putn(best_score); t_puts(": " as *u8)
74 var ok5: i64 = 1
75 if best != 2 { ok5 = 0 }
76 if best_score != 1000 { ok5 = 0 }
77 if ok5 == 1 { pass = pass + 1; t_puts("PASS\n" as *u8) } else { t_puts("FAIL\n" as *u8) }
78
79 // --- T6 CONFIDENT ID requires a margin: the true match is confident, and
80 // a spectrum with no good candidate is NOT ---
81 total = total + 1
82 let conf_true: i64 = msms_confident_id(lib, ang_spec, ang_n, tol)
83 // an empty/garbage observed list: no candidate scores -> not confident
84 let junk: *i64 = sys_mmap(64)
85 junk[0] = 3333333
86 junk[1] = 4444444
87 let conf_junk: i64 = msms_confident_id(lib, junk, 2, tol)
88 t_puts("T6 confident on true spectrum=" as *u8); t_putn(conf_true); t_puts(" on junk spectrum=" as *u8); t_putn(conf_junk); t_puts(" want 1/0: " as *u8)
89 if conf_true == 1 { if conf_junk == 0 { pass = pass + 1; t_puts("PASS\n" as *u8) } else { t_puts("FAIL\n" as *u8) } } else { t_puts("FAIL\n" as *u8) }
90
91 // --- T7 b/y complementarity inside the matcher: every theoretical ion
92 // generated is a real b or y ion of the sequence ---
93 total = total + 1
94 let b3: i64 = pep_b_ion_q4(ang, 3)
95 let y5: i64 = pep_y_ion_q4(ang, 5)
96 let has_b3: i64 = msms_peak_present(b3, ang_spec, ang_n, tol)
97 let has_y5: i64 = msms_peak_present(y5, ang_spec, ang_n, tol)
98 t_puts("T7 generated spectrum contains b3=" as *u8); t_putn(has_b3); t_puts(" y5=" as *u8); t_putn(has_y5); t_puts(" (both present): " as *u8)
99 if has_b3 == 1 { if has_y5 == 1 { pass = pass + 1; t_puts("PASS\n" as *u8) } else { t_puts("FAIL\n" as *u8) } } else { t_puts("FAIL\n" as *u8) }
100
101 // --- T8 a peak NOT in the spectrum is correctly absent (no false match) ---
102 total = total + 1
103 let absent: i64 = msms_peak_present(9999999, ang_spec, ang_n, tol)
104 let present: i64 = msms_peak_present(ang_spec[0], ang_spec, ang_n, tol)
105 t_puts("T8 unrelated m/z present=" as *u8); t_putn(absent); t_puts(" real ion present=" as *u8); t_putn(present); t_puts(" want 0/1: " as *u8)
106 if absent == 0 { if present == 1 { pass = pass + 1; t_puts("PASS\n" as *u8) } else { t_puts("FAIL\n" as *u8) } } else { t_puts("FAIL\n" as *u8) }
107
108 // --- T9 tolerance matters: a peak 0.005 Da off matches at 0.01 tol but
109 // NOT at 0.001 tol -- the matcher respects mass accuracy ---
110 total = total + 1
111 let near: *i64 = sys_mmap(16)
112 near[0] = pep_b_ion_q4(ang, 2) + 50 // +0.005 Da
113 let loose: i64 = msms_peak_present(pep_b_ion_q4(ang, 2), near, 1, 100)
114 let tight: i64 = msms_peak_present(pep_b_ion_q4(ang, 2), near, 1, 10)
115 t_puts("T9 peak +0.005Da: at 0.01 tol=" as *u8); t_putn(loose); t_puts(" at 0.001 tol=" as *u8); t_putn(tight); t_puts(" want 1/0: " as *u8)
116 if loose == 1 { if tight == 0 { pass = pass + 1; t_puts("PASS\n" as *u8) } else { t_puts("FAIL\n" as *u8) } } else { t_puts("FAIL\n" as *u8) }
117
118 // --- T10 fail-closed: a bad sequence refuses spectrum generation +
119 // scores zero ---
120 total = total + 1
121 let bad_spec: *i64 = sys_mmap(64)
122 let bad_n: i64 = msms_theoretical_spectrum("DRVYIHPZ" as *u8, bad_spec, 8)
123 let bad_score: i64 = msms_score_permil("DRVYIHPZ" as *u8, ang_spec, ang_n, tol)
124 let tiny: i64 = msms_theoretical_spectrum("A" as *u8, bad_spec, 8)
125 t_puts("T10 bad-residue spectrum=" as *u8); t_putn(bad_n); t_puts(" score=" as *u8); t_putn(bad_score); t_puts(" 1-mer spectrum=" as *u8); t_putn(tiny); t_puts(" want -1/0/-1: " as *u8)
126 var ok10: i64 = 1
127 if bad_n != MSMS_REFUSED { ok10 = 0 }
128 if bad_score != 0 { ok10 = 0 }
129 if tiny != MSMS_REFUSED { ok10 = 0 }
130 if ok10 == 1 { pass = pass + 1; t_puts("PASS\n" as *u8) } else { t_puts("FAIL\n" as *u8) }
131
132 // --- T11 the SOTA end-to-end: bradykinin's OWN spectrum identifies
133 // bradykinin from the same library, not angiotensin ---
134 total = total + 1
135 let bk_spec: *i64 = sys_mmap(MSMS_MAX_IONS * 8)
136 let bk_n: i64 = msms_theoretical_spectrum(bk, bk_spec, MSMS_MAX_IONS)
137 let bk_best: i64 = msms_best_candidate(lib, bk_spec, bk_n, tol)
138 let bk_conf: i64 = msms_confident_id(lib, bk_spec, bk_n, tol)
139 t_puts("T11 bradykinin spectrum -> candidate " as *u8); t_putn(bk_best); t_puts(" (want 0=bradykinin) confident=" as *u8); t_putn(bk_conf); t_puts(": " as *u8)
140 if bk_best == 0 { if bk_conf == 1 { pass = pass + 1; t_puts("PASS\n" as *u8) } else { t_puts("FAIL\n" as *u8) } } else { t_puts("FAIL\n" as *u8) }
141
142 t_puts("PEPTIDE-MSMS-GATE passed " as *u8); t_putn(pass); t_puts("/" as *u8); t_putn(total)
143 if pass == total { t_puts(" verdict=GREEN\n" as *u8); sys_exit(0); return 0 }
144 t_puts(" verdict=RED\n" as *u8); sys_exit(1); return 1
145}