nx_peptide_identify.nx source
↩ module page · 138 lines · 5574 B
1// nx_peptide_identify.nx -- SOTA CAPSTONE: end-to-end peptide identification
2// from a RAW multiply-charged spectrum. One call chains the whole anchored
3// analytical stack: deconvolve the peak list to a neutral mass, screen that
4// mass against the catalog, and return an identity with a confidence that
5// REFUSES a guess. This is the "spectrum in, identity out" a lab actually
6// wants, and it is the recombination of every rung built before it.
7//
8// THE PIPELINE, EACH STEP ALREADY LIAR-KILLED:
9// 1. DECONVOLVE (nx_peptide_deconv): raw m/z peaks, no known charge -> the
10// neutral monoisotopic mass, plus how many charge-pairs corroborate it.
11// 2. SCREEN (nx_supplement_screen catalog): match that mass against the
12// catalogued peptide actives within a ppm window.
13// 3. JUDGE: an identification is CONFIDENT only if a catalog entry matches
14// AND at least one charge-pair corroborated the deconvolved mass. A
15// single spurious peak, or a mass with no catalog hit, returns
16// "no confident ID" -- it reports the mass it found and declines to name
17// it, rather than forcing a wrong answer.
18//
19// WHY THE REFUSAL MATTERS. A screen that always returns its nearest catalog
20// entry is a liar generator -- it will "identify" noise. The confidence gate
21// here is the difference between an instrument and a rubber stamp: feed it a
22// real BPC-157 ladder and it says BPC-157; feed it noise and it says nothing.
23// The gate proves both.
24//
25// ANCHORED where it is physics: the mass recovery and catalog masses are
26// anchored; the ppm arithmetic is exact. The catalog's regulatory STATUS is
27// asserted (that is nx_supplement_screen's domain), so a caller reads the
28// mass identity as verified and the legal call as advisory -- honestly split.
29//
30// All INTEGER, _q4 = x10^4 Da.
31//
32// Grounding (cited; researcher-groundable):
33// mann_1989_electrospray_charge_deconvolution
34// fda_tainted_supplements_peptide_actives
35// ppm_mass_accuracy_identification_window
36//
37// genealogy_id: peptide_chemistry + analytical_sota + liar_killer
38
39import "nx_syscalls.nx"
40import "nx_peptide.nx"
41import "nx_peptide_deconv.nx"
42import "nx_supplement_screen.nx"
43const IDENT_MAGIC_10000000: i64 = 10000000
44
45// Identification confidence needs BOTH: a catalog hit, and >=1 corroborating
46// charge pair (a real ladder, not one fluke peak).
47const IDENT_MIN_SUPPORT: i64 = 1
48
49// Mass-match window: 50 ppm. At ~1400 Da that is ~0.07 Da = 700 Q4 -- generous
50// enough for the deconvolution's integer rounding, tight enough to reject a
51// different compound (catalog entries differ by many Da).
52const IDENT_PPM_WINDOW: i64 = 50
53
54const IDENT_NO_MATCH: i64 = 0 - 1
55
56struct NxIdentResult {
57 neutral_mass_q4: i64, // deconvolved mass (or DECONV_INVALID)
58 charge_support: i64, // corroborating charge pairs
59 matched_entry: i64, // catalog id or IDENT_NO_MATCH
60 mass_error_ppm_q1: i64, // ppm x10 of the match (or -1)
61 confident: i64, // 1 iff a trustworthy identification
62}
63
64func nx_ident_result_new() -> *NxIdentResult {
65 let r: *NxIdentResult = (sys_mmap(48)) as *NxIdentResult
66 r.neutral_mass_q4 = DECONV_INVALID
67 r.charge_support = 0
68 r.matched_entry = IDENT_NO_MATCH
69 r.mass_error_ppm_q1 = 0 - 1
70 r.confident = 0
71 return r
72}
73
74// ppm x10 between an observed and a reference mass. Q4 scaling cancels.
75func ident_ppm_q1(observed_q4: i64, reference_q4: i64) -> i64 {
76 if reference_q4 <= 0 { return 0 - 1 }
77 var d: i64 = observed_q4 - reference_q4
78 if d < 0 { d = 0 - d }
79 return d * IDENT_MAGIC_10000000 / reference_q4
80}
81
82// Screen a neutral mass against the catalog; return the best (smallest-ppm)
83// entry within the window, or IDENT_NO_MATCH.
84func ident_match_catalog(mass_q4: i64, ppm_window: i64) -> i64 {
85 if mass_q4 <= 0 { return IDENT_NO_MATCH }
86 var i: i64 = 0
87 var best: i64 = IDENT_NO_MATCH
88 var best_ppm: i64 = ppm_window * 10 + 1
89 while i < SUP_N_ENTRY {
90 let em: i64 = sup_mass_q4(i)
91 if em > 0 {
92 let p: i64 = ident_ppm_q1(mass_q4, em)
93 if p >= 0 {
94 if p <= ppm_window * 10 {
95 if p < best_ppm { best_ppm = p; best = i }
96 }
97 }
98 }
99 i = i + 1
100 }
101 return best
102}
103
104// ===== The pipeline ===================================================
105
106func nx_identify_spectrum(peaks: *i64, n: i64, tol_q4: i64) -> *NxIdentResult {
107 let r: *NxIdentResult = nx_ident_result_new()
108 // 1. deconvolve
109 let m: i64 = deconv_from_list(peaks, n, tol_q4)
110 r.neutral_mass_q4 = m
111 if m == DECONV_INVALID { return r }
112 r.charge_support = deconv_support(peaks, n, m, tol_q4)
113 // 2. screen the mass against the catalog
114 let hit: i64 = ident_match_catalog(m, IDENT_PPM_WINDOW)
115 r.matched_entry = hit
116 if hit != IDENT_NO_MATCH {
117 r.mass_error_ppm_q1 = ident_ppm_q1(m, sup_mass_q4(hit))
118 }
119 // 3. judge: confident only with a hit AND corroborating charge support
120 if hit != IDENT_NO_MATCH {
121 if r.charge_support >= IDENT_MIN_SUPPORT { r.confident = 1 }
122 }
123 return r
124}
125
126// Convenience predicates.
127func ident_is_confident(r: *NxIdentResult) -> i64 { return r.confident }
128
129func ident_is_prohibited(r: *NxIdentResult) -> i64 {
130 if r.confident != 1 { return 0 }
131 return sup_is_prohibited(r.matched_entry)
132}
133
134// The identified sequence, or a sentinel if no confident ID.
135func ident_sequence(r: *NxIdentResult) -> *u8 {
136 if r.confident != 1 { return "UNIDENTIFIED" as *u8 }
137 return sup_seq(r.matched_entry)
138}