nx_peptide_isotope.nx source
↩ module page · 135 lines · 5712 B
1// nx_peptide_isotope.nx -- SOTA ANALYTICAL rung: the ISOTOPIC ENVELOPE a
2// high-resolution mass spectrometer actually observes for a peptide, computed
3// from the liar-killed elemental formula (nx_peptide_formula) and NIST natural
4// isotope abundances.
5//
6// WHY THIS IS STATE OF THE ART, NOT A NOMINAL-MASS TOY. A real LC-MS identity
7// confirmation does not match a single mass -- it matches the ISOTOPE PATTERN:
8// the relative heights of the monoisotopic peak (M), M+1, M+2 ... arising from
9// natural 13C/15N/2H/17O/18O/33S/34S. Two compounds can share a nominal mass
10// and be told apart by their envelope. This is the capability that makes the
11// difference between "a peak near 1046" and "yes, that is angiotensin II."
12//
13// THE ENVELOPE IS ANCHORED, NOT ASSERTED. The isotope abundances are NIST
14// measured constants; the pattern is EXACT combinatorics (each element's
15// heavy-isotope contribution, plus the 13C-pair term for M+2). The gate
16// cross-checks the result against an INDEPENDENT empirical law -- "M+1/M is
17// about 1.08% per carbon" -- and against the published ~60% M+1 for a 1 kDa
18// peptide. Computed value == combinatorial value == empirical rule == textbook
19// number is the same triangulation that anchors the mass stack.
20//
21// THE SOTA INSIGHT THIS EXPOSES. For a small molecule the monoisotopic peak
22// is the tallest, so people assume it always is. It is NOT: M+1/M grows ~1.08%
23// per carbon, so past ~90 carbons M+1 exceeds M, and a large peptide's BASE
24// PEAK (the tallest) is M+1 or M+2, not the monoisotopic. iso_base_peak
25// computes which -- a fact a nominal-mass method gets silently wrong.
26//
27// Intensities are relative to the monoisotopic peak = 1000 permil.
28//
29// Grounding (cited; researcher-groundable):
30// nist_atomic_weights_and_isotopic_compositions
31// yergey_1983_polynomial_isotope_distribution
32// published_peptide_isotope_envelopes_1kda
33//
34// genealogy_id: peptide_chemistry + analytical_sota + liar_killer
35
36import "nx_syscalls.nx"
37import "nx_peptide.nx"
38import "nx_peptide_formula.nx"
39const ISO_MAGIC_1000000000: i64 = 1000000000
40
41// ===== Isotope abundance ratios, heavy/light, x10^6 ===================
42// (NIST): the M+1 or M+2 isotope's abundance divided by the most-abundant
43// isotope's abundance, per atom of that element.
44
45const ISO_C13_X6: i64 = 10816 // 13C/12C = 1.07 / 98.93
46const ISO_H2_X6: i64 = 115 // 2H/1H = 0.0115 / 99.9885
47const ISO_N15_X6: i64 = 3653 // 15N/14N = 0.364 / 99.636
48const ISO_O17_X6: i64 = 381 // 17O/16O = 0.038 / 99.757
49const ISO_S33_X6: i64 = 7896 // 33S/32S = 0.75 / 94.99
50
51const ISO_O18_X6: i64 = 2055 // 18O/16O = 0.205 / 99.757 (M+2)
52const ISO_S34_X6: i64 = 44742 // 34S/32S = 4.25 / 94.99 (M+2)
53
54// ===== M+1 relative intensity (permil of the monoisotopic peak) =======
55//
56// Sum over elements of (atom count x heavy/light ratio). Carbon dominates.
57// A ratio x10^6 summed and divided by 1000 gives permil.
58
59func iso_m1_permil(seq: *u8) -> i64 {
60 let f: *NxFormula = pfm_formula(seq)
61 if f.valid != 1 { return PEP_INVALID }
62 var x6: i64 = 0
63 x6 = x6 + f.c * ISO_C13_X6
64 x6 = x6 + f.h * ISO_H2_X6
65 x6 = x6 + f.n * ISO_N15_X6
66 x6 = x6 + f.o * ISO_O17_X6
67 x6 = x6 + f.s * ISO_S33_X6
68 return x6 / 1000
69}
70
71// ===== M+2 relative intensity (permil) ================================
72//
73// Two sources: single +2 isotopes (18O, 34S) AND two independent 13C atoms
74// both being heavy -- the combinatorial term C(nC,2) * (13C/12C)^2, which for
75// a carbon-rich peptide is the LARGEST M+2 contributor. Computed exactly.
76
77func iso_m2_permil(seq: *u8) -> i64 {
78 let f: *NxFormula = pfm_formula(seq)
79 if f.valid != 1 { return PEP_INVALID }
80 // single +2 isotopes, x10^6 -> permil
81 var x6: i64 = 0
82 x6 = x6 + f.o * ISO_O18_X6
83 x6 = x6 + f.s * ISO_S34_X6
84 let single_permil: i64 = x6 / 1000
85 // 13C-pair term: C(nC,2) * (ratio)^2. ratio is x10^6 so ratio^2 is x10^12;
86 // times the pair count, divided by 10^9 gives permil.
87 let nc: i64 = f.c
88 let pairs: i64 = nc * (nc - 1) / 2
89 let c2_permil: i64 = pairs * ISO_C13_X6 * ISO_C13_X6 / ISO_MAGIC_1000000000
90 return single_permil + c2_permil
91}
92
93// ===== The base peak: which is tallest =================================
94//
95// 0 = monoisotopic (M), 1 = M+1, 2 = M+2. For a large peptide this is NOT 0 --
96// the SOTA correction a nominal-mass method misses.
97
98func iso_base_peak(seq: *u8) -> i64 {
99 if pep_seq_valid(seq) != 1 { return PEP_INVALID }
100 let m0: i64 = 1000
101 let m1: i64 = iso_m1_permil(seq)
102 let m2: i64 = iso_m2_permil(seq)
103 var best: i64 = 0
104 var bestv: i64 = m0
105 if m1 > bestv { bestv = m1; best = 1 }
106 if m2 > bestv { bestv = m2; best = 2 }
107 return best
108}
109
110// 1 iff the monoisotopic peak is NOT the tallest -- i.e. a naive "pick the
111// biggest peak = the mass" would report the wrong mass.
112func iso_monoisotopic_is_base_peak(seq: *u8) -> i64 {
113 if iso_base_peak(seq) == 0 { return 1 }
114 return 0
115}
116
117// ===== The independent cross-check (liar-killer) ======================
118//
119// The textbook empirical law: M+1/M (permil) ~= 10.8 per carbon. The full
120// iso_m1_permil must land near this carbon-only estimate, since carbon
121// dominates -- if it did not, the abundance table or the formula is wrong.
122
123func iso_carbon_rule_permil(n_carbon: i64) -> i64 {
124 return n_carbon * ISO_C13_X6 / 1000
125}
126
127// The M+1/M ratio expressed per carbon (permil per carbon), to compare
128// against the ~10.8 empirical constant.
129func iso_m1_per_carbon(seq: *u8) -> i64 {
130 let f: *NxFormula = pfm_formula(seq)
131 if f.valid != 1 { return 0 }
132 if f.c <= 0 { return 0 }
133 let m1: i64 = iso_m1_permil(seq)
134 return m1 * 100 / f.c
135}