nx_chem_isotope_pattern.nx source
↩ module page · 146 lines · 6146 B
1// nx_chem_isotope_pattern.nx -- C2.8b milestone: isotopologue M+1/M+2
2// peak intensity prediction.
3//
4// MS adulterant identification relies on more than exact mass: the
5// natural-abundance isotopologue pattern (M, M+1, M+2 peaks) is a
6// fingerprint of elemental composition.
7//
8// Famous halogen signatures:
9// * Cl in molecule -> M:M+2 ratio approximately 3:1 (75.78% : 24.22%)
10// * Br in molecule -> M:M+2 ratio approximately 1:1 (50.69% : 49.31%)
11// * No Cl/Br -> M+2 small (sub-1% per atom)
12//
13// Sulfur has notable M+2 contribution (~4.5% per S; from 34-S).
14// Carbon has notable M+1 contribution (~1.08% per C; from 13-C).
15//
16// Approach: first-order linear approximation -- each atom contributes
17// its single-substitution probability to M+1 / M+2. Accurate for
18// 0-1 of each halogen + low atom counts; overestimates for multiple
19// halogens (binomial expansion needed; deferred to C2.8b.1).
20//
21// Returns Q4 fraction of M-peak intensity (10000 = 100% of M, so
22// e.g., ethanol M+1 = 232 means M+1 is 2.32% the height of M).
23//
24// Constants pinned (IUPAC 2021 isotopic compositions):
25// ¹H 99.985%, ²H 0.015%
26// ¹²C 98.93%, ¹³C 1.07%
27// ¹⁴N 99.636%, ¹⁵N 0.364%
28// ¹⁶O 99.757%, ¹⁷O 0.038%, ¹⁸O 0.205%
29// ³²S 94.99%, ³³S 0.75%, ³⁴S 4.25%
30// ³⁵Cl 75.78%, ³⁷Cl 24.22%
31// ⁷⁹Br 50.69%, ⁸¹Br 49.31%
32// (¹⁹F, ³¹P, ¹²⁷I monoisotopic -- contribute 0 to M+1/M+2)
33//
34// Honest gaps (deferred):
35// - Multi-halogen binomial expansion: a molecule with 4 Cl atoms
36// gives M+2 = 4 * (37Cl/35Cl) * (35Cl/35Cl)^3 = 4 * 0.3196 = 128%
37// by first-order, but binomial says 4 * 0.3196 * (1 - 0.2422)^3 =
38// ~32% normalized. C2.8b.1 will ship binomial.
39// - M+3, M+4 peaks: relevant for Cl2, Br2, mixed halogen species.
40// C2.8b.2.
41// - Convolution-style isotope-envelope generation a la BRAIN
42// algorithm (Claesen 2012) for spectroscopic-quality envelopes:
43// C2.8c.
44
45import "nx_chem.nx"
46import "nx_chem_molecule.nx"
47
48// Per-atom Q4 fraction-of-M contribution to M+1 peak (parts per 10000).
49const NX_ISO_M1_H: nx_int = 2 // ²H / ¹H = 0.00015 -> 2 in Q4
50const NX_ISO_M1_C: nx_int = 108 // ¹³C / ¹²C = 0.01082 -> 108
51const NX_ISO_M1_N: nx_int = 37 // ¹⁵N / ¹⁴N = 0.00365 -> 37
52const NX_ISO_M1_O: nx_int = 4 // ¹⁷O / ¹⁶O = 0.000381 -> 4
53const NX_ISO_M1_S: nx_int = 79 // ³³S / ³²S = 0.00789 -> 79
54// F (Z=9), P (Z=15), I (Z=53) are monoisotopic; M+1 contribution = 0
55// Cl, Br, Na, K likewise have negligible M+1 (no I+1 isotope or tiny)
56
57// Per-atom Q4 fraction-of-M contribution to M+2 peak.
58const NX_ISO_M2_O: nx_int = 21 // ¹⁸O / ¹⁶O = 0.00206 -> 21
59const NX_ISO_M2_S: nx_int = 447 // ³⁴S / ³²S = 0.0447 -> 447
60const NX_ISO_M2_CL: nx_int = 3196 // ³⁷Cl / ³⁵Cl = 0.3196 -> 3196
61const NX_ISO_M2_BR: nx_int = 9728 // ⁸¹Br / ⁷⁹Br = 0.9728 -> 9728
62// Si has M+2 (3.10%) but rare in pharma; omitted in MVP
63
64// =================================================================
65// Compute Q4 fraction of M-peak intensity that M+1 peak represents.
66// Result = 10000 means M+1 == M; 100 means M+1 is 1% of M.
67//
68// Caller MUST have run nx_chem_compute_implicit_h(m) first so non-
69// bracket atoms have h_count populated (so we count implicit H too).
70//
71// Algorithm: first-order linear sum of per-atom contributions.
72// =================================================================
73func nx_chem_isotope_m1_q4(m: *MolGraph) -> nx_int {
74 var n_c: nx_int = 0
75 var n_h: nx_int = 0 // includes implicit H from h_count
76 var n_n: nx_int = 0
77 var n_o: nx_int = 0
78 var n_s: nx_int = 0
79 var i: nx_int = 0
80 while i < m.n_atoms {
81 let a: *Atom = ((m.atoms as nx_int) + (i * NX_ATOM_BYTES)) as *Atom
82 if a.z == 1 { n_h = n_h + 1 }
83 if a.z == 6 { n_c = n_c + 1 }
84 if a.z == 7 { n_n = n_n + 1 }
85 if a.z == 8 { n_o = n_o + 1 }
86 if a.z == 16 { n_s = n_s + 1 }
87 if a.h_count > 0 { n_h = n_h + a.h_count }
88 i = i + 1
89 }
90 let total: nx_int = (n_c * NX_ISO_M1_C) + (n_h * NX_ISO_M1_H) + (n_n * NX_ISO_M1_N) + (n_o * NX_ISO_M1_O) + (n_s * NX_ISO_M1_S)
91 return total
92}
93
94// =================================================================
95// Compute Q4 fraction of M-peak intensity that M+2 peak represents.
96// =================================================================
97func nx_chem_isotope_m2_q4(m: *MolGraph) -> nx_int {
98 var n_o: nx_int = 0
99 var n_s: nx_int = 0
100 var n_cl: nx_int = 0
101 var n_br: nx_int = 0
102 var i: nx_int = 0
103 while i < m.n_atoms {
104 let a: *Atom = ((m.atoms as nx_int) + (i * NX_ATOM_BYTES)) as *Atom
105 if a.z == 8 { n_o = n_o + 1 }
106 if a.z == 16 { n_s = n_s + 1 }
107 if a.z == 17 { n_cl = n_cl + 1 }
108 if a.z == 35 { n_br = n_br + 1 }
109 i = i + 1
110 }
111 let total: nx_int = (n_o * NX_ISO_M2_O) + (n_s * NX_ISO_M2_S) + (n_cl * NX_ISO_M2_CL) + (n_br * NX_ISO_M2_BR)
112 return total
113}
114
115// =================================================================
116// Halogen signature heuristic: returns
117// 0 if no Cl and no Br
118// 1 if exactly 1 Cl (3:1 M:M+2 signature)
119// 2 if exactly 1 Br (1:1 M:M+2 signature)
120// 3 if mixed / multi-halogen (heavier interpretation needed)
121//
122// Useful for quick adulterant-class detection in MS spectra without
123// requiring fingerprint or substructure match. E.g., sibutramine
124// (1 Cl) gives signature 1; 4-bromosalbutamol (1 Br) gives sig 2.
125// =================================================================
126func nx_chem_halogen_signature(m: *MolGraph) -> nx_int {
127 var n_cl: nx_int = 0
128 var n_br: nx_int = 0
129 var i: nx_int = 0
130 while i < m.n_atoms {
131 let a: *Atom = ((m.atoms as nx_int) + (i * NX_ATOM_BYTES)) as *Atom
132 if a.z == 17 { n_cl = n_cl + 1 }
133 if a.z == 35 { n_br = n_br + 1 }
134 i = i + 1
135 }
136 if n_cl == 0 {
137 if n_br == 0 { return 0 }
138 if n_br == 1 { return 2 }
139 return 3
140 }
141 if n_cl == 1 {
142 if n_br == 0 { return 1 }
143 return 3
144 }
145 return 3
146}