code wiki / (root) / nx_peptide_identify_test.nx

nx_peptide_identify_test.nx source

↩ module page · 156 lines · 8945 B

1// nx_peptide_identify_test.nx -- gate for the end-to-end identification 2// capstone. The liar-kill is the pair T1/T4: a REAL spectrum identifies the 3// right prohibited peptide confidently, and NOISE identifies nothing. A screen 4// that named noise would be a liar generator; this refuses. 5// expect_exit: 0 license_tier: ORIGINAL 6import "nx_syscalls.nx" 7import "nx_peptide.nx" 8import "nx_peptide_deconv.nx" 9import "nx_supplement_screen.nx" 10import "nx_peptide_identify.nx" 11 12func 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 } 13func 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 } 14 15func t_streq(a: *u8, b: *u8) -> i64 { var i: i64 = 0; while 1 == 1 { if a[i] != b[i] { return 0 } if a[i] == (0 as u8) { return 1 } i = i + 1 } return 1 } 16 17// Synthesize a real charge ladder (z=2,3,4) for a catalog entry, ascending. 18func make_ladder(id: i64, out: *i64) -> i64 { 19 out[0] = sup_mz_q4(id, 4) 20 out[1] = sup_mz_q4(id, 3) 21 out[2] = sup_mz_q4(id, 2) 22 return 3 23} 24 25func main() -> i64 { 26 var pass: i64 = 0 27 var total: i64 = 0 28 let tol: i64 = 100 29 30 // --- T1 THE HEADLINE: a real BPC-157 spectrum -> confident BPC-157 ID, 31 // flagged prohibited, low ppm error --- 32 total = total + 1 33 let bpc_spec: *i64 = sys_mmap(64) 34 let bn: i64 = make_ladder(SUP_BPC157, bpc_spec) 35 let r1: *NxIdentResult = nx_identify_spectrum(bpc_spec, bn, tol) 36 t_puts("T1 BPC-157 spectrum -> entry=" as *u8); t_putn(r1.matched_entry); t_puts(" mass=" as *u8); t_putn(r1.neutral_mass_q4); t_puts(" support=" as *u8); t_putn(r1.charge_support); t_puts(" ppm_x10=" as *u8); t_putn(r1.mass_error_ppm_q1); t_puts(" confident=" as *u8); t_putn(r1.confident); t_puts(" prohibited=" as *u8); t_putn(ident_is_prohibited(r1)); t_puts(": " as *u8) 37 var ok1: i64 = 1 38 if r1.matched_entry != SUP_BPC157 { ok1 = 0 } 39 if r1.confident != 1 { ok1 = 0 } 40 if ident_is_prohibited(r1) != 1 { ok1 = 0 } 41 if r1.mass_error_ppm_q1 < 0 { ok1 = 0 } 42 if r1.mass_error_ppm_q1 > 300 { ok1 = 0 } // within 30 ppm 43 if ok1 == 1 { pass = pass + 1; t_puts("PASS\n" as *u8) } else { t_puts("FAIL\n" as *u8) } 44 45 // --- T2 the identified SEQUENCE is BPC-157's, read back from the ID --- 46 total = total + 1 47 let seq: *u8 = ident_sequence(r1) 48 t_puts("T2 identified sequence = " as *u8); t_puts(seq); t_puts(" want GEPPPGKPADDAGLV: " as *u8) 49 if t_streq(seq, "GEPPPGKPADDAGLV" as *u8) == 1 { pass = pass + 1; t_puts("PASS\n" as *u8) } else { t_puts("FAIL\n" as *u8) } 50 51 // --- T3 a SECOND compound (angiotensin) identifies correctly, so the 52 // pipeline is not hard-wired to one answer --- 53 total = total + 1 54 let ang_spec: *i64 = sys_mmap(64) 55 let an: i64 = make_ladder(SUP_ANGII, ang_spec) 56 let r3: *NxIdentResult = nx_identify_spectrum(ang_spec, an, tol) 57 t_puts("T3 angiotensin spectrum -> entry=" as *u8); t_putn(r3.matched_entry); t_puts(" (want " as *u8); t_putn(SUP_ANGII); t_puts(") confident=" as *u8); t_putn(r3.confident); t_puts(": " as *u8) 58 var ok3: i64 = 1 59 if r3.matched_entry != SUP_ANGII { ok3 = 0 } 60 if r3.confident != 1 { ok3 = 0 } 61 if ok3 == 1 { pass = pass + 1; t_puts("PASS\n" as *u8) } else { t_puts("FAIL\n" as *u8) } 62 63 // --- T4 THE LIAR-KILL: NOISE identifies NOTHING. Random peaks that are 64 // not a real charge ladder -> no confident ID --- 65 total = total + 1 66 let noise: *i64 = sys_mmap(64) 67 noise[0] = 3210000 68 noise[1] = 5550000 69 noise[2] = 8880000 70 let r4: *NxIdentResult = nx_identify_spectrum(noise, 3, tol) 71 t_puts("T4 NOISE spectrum -> confident=" as *u8); t_putn(r4.confident); t_puts(" entry=" as *u8); t_putn(r4.matched_entry); t_puts(" (want 0/-1: names nothing): " as *u8) 72 var ok4: i64 = 1 73 if r4.confident != 0 { ok4 = 0 } 74 if ok4 == 1 { pass = pass + 1; t_puts("PASS\n" as *u8) } else { t_puts("FAIL\n" as *u8) } 75 76 // --- T5 a peptide NOT in the catalog: deconvolves to a MASS but declines 77 // to NAME it (honest: found a mass, no ID) --- 78 total = total + 1 79 // GHRP-2 backbone-ish not in catalog; use a random real peptide's ladder. 80 let unk: *i64 = sys_mmap(64) 81 let unkM: i64 = pep_mass_mono_q4("KKKKKKKKKK" as *u8) // polylysine, not catalogued 82 unk[0] = pep_mz_q4(unkM, 3) 83 unk[1] = pep_mz_q4(unkM, 2) 84 unk[2] = pep_mz_q4(unkM, 1) 85 let r5: *NxIdentResult = nx_identify_spectrum(unk, 3, tol) 86 t_puts("T5 uncatalogued peptide -> mass found=" as *u8); t_putn(r5.neutral_mass_q4); t_puts(" (>0) entry=" as *u8); t_putn(r5.matched_entry); t_puts(" confident=" as *u8); t_putn(r5.confident); t_puts(" (mass yes, ID no): " as *u8) 87 var ok5: i64 = 1 88 if r5.neutral_mass_q4 <= 0 { ok5 = 0 } // it DID recover a mass 89 if r5.matched_entry != IDENT_NO_MATCH { ok5 = 0 } // but no catalog name 90 if r5.confident != 0 { ok5 = 0 } 91 if ok5 == 1 { pass = pass + 1; t_puts("PASS\n" as *u8) } else { t_puts("FAIL\n" as *u8) } 92 93 // --- T6 a SINGLE peak cannot be deconvolved -> no mass, no ID --- 94 total = total + 1 95 let onep: *i64 = sys_mmap(16) 96 onep[0] = sup_mz_q4(SUP_BPC157, 2) 97 let r6: *NxIdentResult = nx_identify_spectrum(onep, 1, tol) 98 t_puts("T6 single peak -> mass=" as *u8); t_putn(r6.neutral_mass_q4); t_puts(" confident=" as *u8); t_putn(r6.confident); t_puts(" want -1/0: " as *u8) 99 var ok6: i64 = 1 100 if r6.neutral_mass_q4 != DECONV_INVALID { ok6 = 0 } 101 if r6.confident != 0 { ok6 = 0 } 102 if ok6 == 1 { pass = pass + 1; t_puts("PASS\n" as *u8) } else { t_puts("FAIL\n" as *u8) } 103 104 // --- T7 ppm arithmetic sanity + catalog matcher direct --- 105 total = total + 1 106 let exact: i64 = ident_ppm_q1(14187044, 14187044) 107 let off: i64 = ident_ppm_q1(14187044 + 700, 14187044) // +0.07 Da ~50 ppm 108 let match_bpc: i64 = ident_match_catalog(14187044, IDENT_PPM_WINDOW) 109 let no_match: i64 = ident_match_catalog(50000000, IDENT_PPM_WINDOW) // 5000 Da, nothing 110 t_puts("T7 ppm exact=" as *u8); t_putn(exact); t_puts(" +0.07Da=" as *u8); t_putn(off); t_puts(" x10 catalog-match(BPC mass)=" as *u8); t_putn(match_bpc); t_puts(" match(5000Da)=" as *u8); t_putn(no_match); t_puts(": " as *u8) 111 var ok7: i64 = 1 112 if exact != 0 { ok7 = 0 } 113 if off < 400 { ok7 = 0 } 114 if off > 550 { ok7 = 0 } 115 if match_bpc != SUP_BPC157 { ok7 = 0 } 116 if no_match != IDENT_NO_MATCH { ok7 = 0 } 117 if ok7 == 1 { pass = pass + 1; t_puts("PASS\n" as *u8) } else { t_puts("FAIL\n" as *u8) } 118 119 // --- T8 charge support: a 3-peak real ladder corroborates with 2 pairs; 120 // a 2-peak ladder with 1 -- more charge states = more confidence --- 121 total = total + 1 122 let two: *i64 = sys_mmap(32) 123 two[0] = sup_mz_q4(SUP_BPC157, 3) 124 two[1] = sup_mz_q4(SUP_BPC157, 2) 125 let r8: *NxIdentResult = nx_identify_spectrum(two, 2, tol) 126 t_puts("T8 2-peak ladder support=" as *u8); t_putn(r8.charge_support); t_puts(" (1) vs 3-peak support=" as *u8); t_putn(r1.charge_support); t_puts(" (2), both confident: " as *u8) 127 var ok8: i64 = 1 128 if r8.charge_support != 1 { ok8 = 0 } 129 if r1.charge_support != 2 { ok8 = 0 } 130 if r8.confident != 1 { ok8 = 0 } 131 if ok8 == 1 { pass = pass + 1; t_puts("PASS\n" as *u8) } else { t_puts("FAIL\n" as *u8) } 132 133 // --- T9 TB-500 and GHRP-6 also identify from their own spectra --- 134 total = total + 1 135 let tb: *i64 = sys_mmap(64) 136 make_ladder(SUP_TB500, tb) 137 let r_tb: *NxIdentResult = nx_identify_spectrum(tb, 3, tol) 138 let gh: *i64 = sys_mmap(64) 139 make_ladder(SUP_GHRP6, gh) 140 let r_gh: *NxIdentResult = nx_identify_spectrum(gh, 3, tol) 141 t_puts("T9 TB-500 spectrum -> entry=" as *u8); t_putn(r_tb.matched_entry); t_puts(" GHRP-6 -> entry=" as *u8); t_putn(r_gh.matched_entry); t_puts(" want 1/2: " as *u8) 142 var ok9: i64 = 1 143 if r_tb.matched_entry != SUP_TB500 { ok9 = 0 } 144 if r_gh.matched_entry != SUP_GHRP6 { ok9 = 0 } 145 if ok9 == 1 { pass = pass + 1; t_puts("PASS\n" as *u8) } else { t_puts("FAIL\n" as *u8) } 146 147 // --- T10 UNIDENTIFIED sequence sentinel when no confident ID --- 148 total = total + 1 149 let noise_seq: *u8 = ident_sequence(r4) 150 t_puts("T10 noise ID sequence = " as *u8); t_puts(noise_seq); t_puts(" want UNIDENTIFIED: " as *u8) 151 if t_streq(noise_seq, "UNIDENTIFIED" as *u8) == 1 { pass = pass + 1; t_puts("PASS\n" as *u8) } else { t_puts("FAIL\n" as *u8) } 152 153 t_puts("PEPTIDE-IDENTIFY-GATE passed " as *u8); t_putn(pass); t_puts("/" as *u8); t_putn(total) 154 if pass == total { t_puts(" verdict=GREEN\n" as *u8); sys_exit(0); return 0 } 155 t_puts(" verdict=RED\n" as *u8); sys_exit(1); return 1 156}