nx_peptide_formula.nx source
↩ module page · 176 lines · 7768 B
1// nx_peptide_formula.nx -- LIAR-KILLER for the peptide mass stack. Computes a
2// peptide's monoisotopic mass by a SECOND, INDEPENDENT path: counting the
3// C/H/N/O/S atoms of every residue and multiplying by CODATA element masses,
4// rather than summing a pre-rounded residue-mass table (nx_peptide).
5//
6// WHY THIS IS REAL LIAR-KILLING, NOT SELF-CONSISTENCY. nx_peptide sums a
7// table of residue masses. If that table had a typo, its gate -- which checks
8// the sum against a literature number -- would still catch a SINGLE-peptide
9// error, but a systematic table error could hide. This organ derives the same
10// masses from a DIFFERENT primitive (atom counts x element masses), and the
11// gate demands three things agree:
12// (1) this elemental mass == nx_peptide's residue-sum mass (two methods)
13// (2) both == the PUBLISHED monoisotopic mass (external)
14// (3) the derived MOLECULAR FORMULA == the PUBLISHED formula (external)
15// For a table typo to survive, the SAME error would have to appear in the
16// residue-mass table AND the atom-count table AND match the published formula
17// AND the published mass -- which is not a mistake, it is a conspiracy. Two
18// independent representations converging on external ground truth is the
19// definition of a killed liar.
20//
21// The molecular-formula check (3) is the sharpest: angiotensin II is the
22// published C50H71N13O12, bradykinin C50H73N15O11. Those are hand-checkable
23// integers, not a mass that hides rounding. If the atom table is wrong the
24// formula is visibly wrong.
25//
26// All INTEGER, _q4 = x10^4 Da, matching nx_peptide.
27//
28// Grounding (cited; researcher-groundable):
29// codata_monoisotopic_atomic_masses_H_C_N_O_S
30// published_molecular_formulas_angiotensin_ii_bradykinin
31// iupac_residue_elemental_compositions
32//
33// genealogy_id: peptide_chemistry + liar_killer
34
35import "nx_syscalls.nx"
36import "nx_peptide.nx"
37
38// ===== Element monoisotopic masses, Q6 (x10^6 Da) =====================
39// Most-abundant isotope, CODATA/AME2020. C is the 12C anchor = 12.000000.
40// HELD AT Q6, NOT Q4: per-atom Q4 rounding accumulates ~0.1 Q4 per atom, so a
41// 146-atom peptide drifts ~16 Q4 low -- which is exactly the discrepancy the
42// liar-killer gate exposed against the residue-sum method. Summing at Q6 and
43// rounding to Q4 ONCE at the end makes the two independent methods agree.
44
45const PFM_H_Q6: i64 = 1007825 // 1.00782503
46const PFM_C_Q6: i64 = 12000000 // 12.00000000
47const PFM_N_Q6: i64 = 14003074 // 14.00307401
48const PFM_O_Q6: i64 = 15994915 // 15.99491462
49const PFM_S_Q6: i64 = 31972071 // 31.97207069
50
51// Water = 2H + O, in Q4 for the residue-stripping check.
52const PFM_WATER_Q4: i64 = 180106 // (2*1007825 + 15994915 + 50)/100
53
54// ===== Per-residue elemental composition (residue = amino acid - H2O) ==
55// Returned packed so one lookup yields all five counts without 5 branches.
56// Encoding: C*1 + H*100 + N*10000 + O*1000000 + S*100000000. Each field <100.
57
58func pfm_atoms_packed(aa: i64) -> i64 {
59 if aa == PEP_G { return 2 + 3*100 + 1*10000 + 1*1000000 + 0 } // C2H3NO
60 if aa == PEP_A { return 3 + 5*100 + 1*10000 + 1*1000000 + 0 } // C3H5NO
61 if aa == PEP_S { return 3 + 5*100 + 1*10000 + 2*1000000 + 0 } // C3H5NO2
62 if aa == PEP_P { return 5 + 7*100 + 1*10000 + 1*1000000 + 0 } // C5H7NO
63 if aa == PEP_V { return 5 + 9*100 + 1*10000 + 1*1000000 + 0 } // C5H9NO
64 if aa == PEP_T { return 4 + 7*100 + 1*10000 + 2*1000000 + 0 } // C4H7NO2
65 if aa == PEP_C { return 3 + 5*100 + 1*10000 + 1*1000000 + 1*100000000 }// C3H5NOS
66 if aa == PEP_L { return 6 + 11*100 + 1*10000 + 1*1000000 + 0 } // C6H11NO
67 if aa == PEP_I { return 6 + 11*100 + 1*10000 + 1*1000000 + 0 } // C6H11NO
68 if aa == PEP_N { return 4 + 6*100 + 2*10000 + 2*1000000 + 0 } // C4H6N2O2
69 if aa == PEP_D { return 4 + 5*100 + 1*10000 + 3*1000000 + 0 } // C4H5NO3
70 if aa == PEP_Q { return 5 + 8*100 + 2*10000 + 2*1000000 + 0 } // C5H8N2O2
71 if aa == PEP_K { return 6 + 12*100 + 2*10000 + 1*1000000 + 0 } // C6H12N2O
72 if aa == PEP_E { return 5 + 7*100 + 1*10000 + 3*1000000 + 0 } // C5H7NO3
73 if aa == PEP_M { return 5 + 9*100 + 1*10000 + 1*1000000 + 1*100000000 }// C5H9NOS
74 if aa == PEP_H { return 6 + 7*100 + 3*10000 + 1*1000000 + 0 } // C6H7N3O
75 if aa == PEP_F { return 9 + 9*100 + 1*10000 + 1*1000000 + 0 } // C9H9NO
76 if aa == PEP_R { return 6 + 12*100 + 4*10000 + 1*1000000 + 0 } // C6H12N4O
77 if aa == PEP_Y { return 9 + 9*100 + 1*10000 + 2*1000000 + 0 } // C9H9NO2
78 if aa == PEP_W { return 11 + 10*100 + 2*10000 + 1*1000000 + 0 } // C11H10N2O
79 return 0 - 1
80}
81
82func pfm_unpack_c(p: i64) -> i64 { return p % 100 }
83func pfm_unpack_h(p: i64) -> i64 { return (p / 100) % 100 }
84func pfm_unpack_n(p: i64) -> i64 { return (p / 10000) % 100 }
85func pfm_unpack_o(p: i64) -> i64 { return (p / 1000000) % 100 }
86func pfm_unpack_s(p: i64) -> i64 { return (p / 100000000) % 100 }
87
88// ===== Molecular formula of a whole peptide ===========================
89//
90// Sum residue atoms + one water (adds 2 H and 1 O for the chain termini).
91// Returns the packed total; unpack for each element. REFUSES on a bad residue.
92
93struct NxFormula {
94 c: i64,
95 h: i64,
96 n: i64,
97 o: i64,
98 s: i64,
99 valid: i64,
100}
101
102func pfm_formula(seq: *u8) -> *NxFormula {
103 let f: *NxFormula = (sys_mmap(48)) as *NxFormula
104 f.c = 0
105 f.h = 0
106 f.n = 0
107 f.o = 0
108 f.s = 0
109 f.valid = 1
110 if pep_seq_valid(seq) != 1 { f.valid = 0; return f }
111 let n: i64 = pep_len(seq)
112 var i: i64 = 0
113 while i < n {
114 let c: i64 = seq[i] as i64
115 let aa: i64 = pep_aa_from_char(c)
116 let p: i64 = pfm_atoms_packed(aa)
117 if p < 0 { f.valid = 0 }
118 if p >= 0 {
119 f.c = f.c + pfm_unpack_c(p)
120 f.h = f.h + pfm_unpack_h(p)
121 f.n = f.n + pfm_unpack_n(p)
122 f.o = f.o + pfm_unpack_o(p)
123 f.s = f.s + pfm_unpack_s(p)
124 }
125 i = i + 1
126 }
127 // + water: 2 H, 1 O.
128 f.h = f.h + 2
129 f.o = f.o + 1
130 return f
131}
132
133// ===== Mass from the elemental formula (the independent path) ==========
134
135func pfm_mass_q4(seq: *u8) -> i64 {
136 let f: *NxFormula = pfm_formula(seq)
137 if f.valid != 1 { return PEP_INVALID }
138 var q6: i64 = 0
139 q6 = q6 + f.c * PFM_C_Q6
140 q6 = q6 + f.h * PFM_H_Q6
141 q6 = q6 + f.n * PFM_N_Q6
142 q6 = q6 + f.o * PFM_O_Q6
143 q6 = q6 + f.s * PFM_S_Q6
144 // Round Q6 -> Q4 ONCE, at the end, so per-atom rounding cannot accumulate.
145 return (q6 + 50) / 100
146}
147
148// ===== The cross-check ================================================
149//
150// 1 iff the elemental-path mass and nx_peptide's residue-sum mass agree
151// within tol_q4 -- two independent methods on the same molecule. A rounding
152// budget of a few Q4 units is expected (the residue table and the element
153// table round independently); a real discrepancy is a table bug.
154func pfm_methods_agree(seq: *u8, tol_q4: i64) -> i64 {
155 let a: i64 = pfm_mass_q4(seq)
156 let b: i64 = pep_mass_mono_q4(seq)
157 if a == PEP_INVALID { return 0 }
158 if b == PEP_INVALID { return 0 }
159 var d: i64 = a - b
160 if d < 0 { d = 0 - d }
161 if d <= tol_q4 { return 1 }
162 return 0
163}
164
165// 1 iff the derived formula equals a published (c,h,n,o,s) tuple -- the
166// sharpest external check, on hand-verifiable integers.
167func pfm_formula_matches(seq: *u8, c: i64, h: i64, n: i64, o: i64, s: i64) -> i64 {
168 let f: *NxFormula = pfm_formula(seq)
169 if f.valid != 1 { return 0 }
170 if f.c != c { return 0 }
171 if f.h != h { return 0 }
172 if f.n != n { return 0 }
173 if f.o != o { return 0 }
174 if f.s != s { return 0 }
175 return 1
176}