code wiki / (root) / nx_peptide_denovo_test.nx

nx_peptide_denovo_test.nx source

↩ module page · 162 lines · 8885 B

1// nx_peptide_denovo_test.nx -- gate for de novo sequencing. The liar-kill is 2// the ROUND TRIP (T1): take a known peptide, generate its real b-ion ladder, 3// de novo sequence it back with NO knowledge of the answer, and recover the 4// original sequence. If any residue-mass gap were wrong, a letter would be 5// wrong. Bradykinin (no Leu/Ile) gives an EXACT recovery; angiotensin 6// documents the honest Leu/Ile ambiguity. 7// expect_exit: 0 license_tier: ORIGINAL 8import "nx_syscalls.nx" 9import "nx_peptide.nx" 10import "nx_peptide_denovo.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// Build the full b-ion ladder (b1..b(n-1)) of a sequence into out. Returns n-1. 18func build_b_ladder(seq: *u8, out: *i64) -> i64 { 19 let n: i64 = pep_len(seq) 20 var i: i64 = 1 21 while i < n { 22 out[i - 1] = pep_b_ion_q4(seq, i) 23 i = i + 1 24 } 25 return n - 1 26} 27 28func main() -> i64 { 29 var pass: i64 = 0 30 var total: i64 = 0 31 let tol: i64 = 50 // 0.005 Da 32 33 // --- T1 THE LIAR-KILL: round-trip bradykinin (no Leu/Ile) EXACTLY --- 34 total = total + 1 35 let bk: *u8 = "RPPGFSPFR" as *u8 36 let bk_M: i64 = pep_mass_mono_q4(bk) 37 let bk_ladder: *i64 = sys_mmap(128 * 8) 38 let bk_nb: i64 = build_b_ladder(bk, bk_ladder) 39 let out1: *u8 = sys_mmap(128) 40 let len1: i64 = denovo_sequence(bk_ladder, bk_nb, bk_M, tol, out1) 41 t_puts("T1 de novo bradykinin ladder -> '" as *u8); t_puts(out1); t_puts("' (len " as *u8); t_putn(len1); t_puts(") vs RPPGFSPFR: " as *u8) 42 var ok1: i64 = 1 43 if len1 != 9 { ok1 = 0 } 44 if t_streq(out1, "RPPGFSPFR" as *u8) != 1 { ok1 = 0 } 45 if ok1 == 1 { pass = pass + 1; t_puts("PASS\n" as *u8) } else { t_puts("FAIL\n" as *u8) } 46 47 // --- T2 a SECOND peptide round-trips (not hard-coded) --- 48 total = total + 1 49 let gh: *u8 = "HWAWFK" as *u8 // GHRP-6 backbone, no Leu/Ile 50 let gh_M: i64 = pep_mass_mono_q4(gh) 51 let gh_ladder: *i64 = sys_mmap(128 * 8) 52 let gh_nb: i64 = build_b_ladder(gh, gh_ladder) 53 let out2: *u8 = sys_mmap(128) 54 denovo_sequence(gh_ladder, gh_nb, gh_M, tol, out2) 55 t_puts("T2 de novo '" as *u8); t_puts(out2); t_puts("' vs HWAWFK: " as *u8) 56 if t_streq(out2, "HWAWFK" as *u8) == 1 { pass = pass + 1; t_puts("PASS\n" as *u8) } else { t_puts("FAIL\n" as *u8) } 57 58 // --- T3 THE HONEST AMBIGUITY: angiotensin has an Ile; de novo reads it 59 // as canonical L and FLAGS the position -- not a silent wrong call --- 60 total = total + 1 61 let ang: *u8 = "DRVYIHPF" as *u8 62 let ang_M: i64 = pep_mass_mono_q4(ang) 63 let ang_ladder: *i64 = sys_mmap(128 * 8) 64 let ang_nb: i64 = build_b_ladder(ang, ang_ladder) 65 let out3: *u8 = sys_mmap(128) 66 denovo_sequence(ang_ladder, ang_nb, ang_M, tol, out3) 67 let amb: i64 = denovo_ambiguous_count(ang_ladder, ang_nb, ang_M, tol) 68 t_puts("T3 de novo angiotensin -> '" as *u8); t_puts(out3); t_puts("' (Ile read as L) ambiguous positions=" as *u8); t_putn(amb); t_puts(" want 'DRVYLHPF'/1: " as *u8) 69 var ok3: i64 = 1 70 if t_streq(out3, "DRVYLHPF" as *u8) != 1 { ok3 = 0 } // I recovered as canonical L 71 if amb != 1 { ok3 = 0 } // exactly one L/I ambiguity 72 if ok3 == 1 { pass = pass + 1; t_puts("PASS\n" as *u8) } else { t_puts("FAIL\n" as *u8) } 73 74 // --- T4 residue-from-gap is exact for a battery of residues --- 75 total = total + 1 76 let rD: i64 = denovo_residue(pep_residue_mono_q4(PEP_D), tol) 77 let rW: i64 = denovo_residue(pep_residue_mono_q4(PEP_W), tol) 78 let rG: i64 = denovo_residue(pep_residue_mono_q4(PEP_G), tol) 79 let rR: i64 = denovo_residue(pep_residue_mono_q4(PEP_R), tol) 80 t_puts("T4 residue-from-mass D/W/G/R = " as *u8); t_putn(rD); t_puts("/" as *u8); t_putn(rW); t_puts("/" as *u8); t_putn(rG); t_puts("/" as *u8); t_putn(rR); t_puts(" want 2/18/5/14: " as *u8) 81 var ok4: i64 = 1 82 if rD != PEP_D { ok4 = 0 } 83 if rW != PEP_W { ok4 = 0 } 84 if rG != PEP_G { ok4 = 0 } 85 if rR != PEP_R { ok4 = 0 } 86 if ok4 == 1 { pass = pass + 1; t_puts("PASS\n" as *u8) } else { t_puts("FAIL\n" as *u8) } 87 88 // --- T5 Leu/Ile isobaric: both masses resolve to the SAME canonical L, 89 // and the position is flagged ambiguous --- 90 total = total + 1 91 let mL: i64 = pep_residue_mono_q4(PEP_L) 92 let mI: i64 = pep_residue_mono_q4(PEP_I) 93 let rL: i64 = denovo_residue(mL, tol) 94 let rI: i64 = denovo_residue(mI, tol) 95 let flagL: i64 = denovo_is_leu_ile(mL, tol) 96 let flagD: i64 = denovo_is_leu_ile(pep_residue_mono_q4(PEP_D), tol) 97 t_puts("T5 Leu mass->" as *u8); t_putn(rL); t_puts(" Ile mass->" as *u8); t_putn(rI); t_puts(" (both canonical L=9) flag L/I=" as *u8); t_putn(flagL); t_puts(" flag D=" as *u8); t_putn(flagD); t_puts(": " as *u8) 98 var ok5: i64 = 1 99 if mL != mI { ok5 = 0 } // truly isobaric 100 if rL != PEP_L { ok5 = 0 } 101 if rI != PEP_L { ok5 = 0 } // Ile also resolves to canonical L 102 if flagL != 1 { ok5 = 0 } 103 if flagD != 0 { ok5 = 0 } 104 if ok5 == 1 { pass = pass + 1; t_puts("PASS\n" as *u8) } else { t_puts("FAIL\n" as *u8) } 105 106 // --- T6 denovo_char is the exact inverse of pep_aa_from_char --- 107 total = total + 1 108 var aa6: i64 = 0 109 var ok6: i64 = 1 110 while aa6 < PEP_N_AA { 111 let ch: i64 = denovo_char(aa6) 112 let back: i64 = pep_aa_from_char(ch) 113 if back != aa6 { ok6 = 0 } 114 aa6 = aa6 + 1 115 } 116 t_puts("T6 denovo_char inverts pep_aa_from_char for all 20 residues: " as *u8) 117 if ok6 == 1 { pass = pass + 1; t_puts("PASS\n" as *u8) } else { t_puts("FAIL\n" as *u8) } 118 119 // --- T7 a gap matching NO residue -> UNKNOWN (fail-closed, no false call) --- 120 total = total + 1 121 let junk_gap: i64 = denovo_residue(500000, tol) // 50 Da: no residue 122 let big_gap: i64 = denovo_residue(9999999, tol) 123 t_puts("T7 unmatchable gaps 50Da/1000Da = " as *u8); t_putn(junk_gap); t_puts("/" as *u8); t_putn(big_gap); t_puts(" want -1/-1: " as *u8) 124 var ok7: i64 = 1 125 if junk_gap != DENOVO_UNKNOWN { ok7 = 0 } 126 if big_gap != DENOVO_UNKNOWN { ok7 = 0 } 127 if ok7 == 1 { pass = pass + 1; t_puts("PASS\n" as *u8) } else { t_puts("FAIL\n" as *u8) } 128 129 // --- T8 a corrupt ladder (one ion shifted off any residue) REFUSES to 130 // sequence -- a wrong spectrum yields no answer, not a wrong one --- 131 total = total + 1 132 let corrupt: *i64 = sys_mmap(128 * 8) 133 var c8: i64 = 0 134 while c8 < bk_nb { corrupt[c8] = bk_ladder[c8]; c8 = c8 + 1 } 135 corrupt[2] = corrupt[2] + 500000 // shift b3 by 50 Da: gap no longer a residue 136 let outc: *u8 = sys_mmap(128) 137 let lenc: i64 = denovo_sequence(corrupt, bk_nb, bk_M, tol, outc) 138 t_puts("T8 corrupt ladder -> len=" as *u8); t_putn(lenc); t_puts(" want -1 (refuses): " as *u8) 139 if lenc == DENOVO_UNKNOWN { pass = pass + 1; t_puts("PASS\n" as *u8) } else { t_puts("FAIL\n" as *u8) } 140 141 // --- T9 the reconstructed sequence, fed BACK, regenerates the same mass 142 // (self-consistent round trip of the round trip) --- 143 total = total + 1 144 let rt_mass: i64 = pep_mass_mono_q4(out1) // out1 = "RPPGFSPFR" 145 t_puts("T9 de novo output mass=" as *u8); t_putn(rt_mass); t_puts(" vs original bradykinin " as *u8); t_putn(bk_M); t_puts(": " as *u8) 146 if rt_mass == bk_M { pass = pass + 1; t_puts("PASS\n" as *u8) } else { t_puts("FAIL\n" as *u8) } 147 148 // --- T10 tolerance matters: at ULTRA-tight tol a real gap still matches, 149 // at absurdly tight tol (1 Q4) rounding blocks it (fail-closed) --- 150 total = total + 1 151 let ok_tol: i64 = denovo_residue(pep_residue_mono_q4(PEP_F), 50) 152 let tight_tol: i64 = denovo_residue(pep_residue_mono_q4(PEP_F) + 30, 5) // +0.003 Da, 0.0005 tol 153 t_puts("T10 Phe at 0.005 tol=" as *u8); t_putn(ok_tol); t_puts(" (=4) Phe+0.003Da at 0.0005 tol=" as *u8); t_putn(tight_tol); t_puts(" (=-1): " as *u8) 154 var ok10: i64 = 1 155 if ok_tol != PEP_F { ok10 = 0 } 156 if tight_tol != DENOVO_UNKNOWN { ok10 = 0 } 157 if ok10 == 1 { pass = pass + 1; t_puts("PASS\n" as *u8) } else { t_puts("FAIL\n" as *u8) } 158 159 t_puts("PEPTIDE-DENOVO-GATE passed " as *u8); t_putn(pass); t_puts("/" as *u8); t_putn(total) 160 if pass == total { t_puts(" verdict=GREEN\n" as *u8); sys_exit(0); return 0 } 161 t_puts(" verdict=RED\n" as *u8); sys_exit(1); return 1 162}