nx_peptide.nx source

1// nx_peptide.nx -- CHEMISTRY SUITE / PEPTIDE PRIMARY-STRUCTURE rung. Turns 2// an amino-acid SEQUENCE into real chemistry: exact mass, charge state, 3// isoelectric point, UV absorbance, hydropathy, and the b/y fragment-ion 4// ladder an MS/MS run actually observes. 5// 6// This is the missing rung between the ecosystem's molecule stack 7// (nx_chem_molecule / nx_chem_mass / nx_chem_isotope_pattern / nx_chem_peak_list) 8// and peptide work: the mass stack could weigh a MolGraph but nothing could 9// go from "DRVYIHPF" to 1045.5345 Da, so the NIST Peptide Tandem Mass 10// Spectral Libraries (SRD 1c) in the corpus had no local counterpart to 11// match against. This closes that. 12// 13// ALL INTEGER, no float: 14// *_q4 = value x 10^4 (mass in Da -> micro-AMU; matches nx_chem_mass) 15// *_q3 = value x 10^3 (charge, hydropathy, decade-log intermediates) 16// *_milli = value x 10^3 (pH and pKa) 17// 18// THE exceed: a sequence viewer prints letters. This computes the mass an 19// instrument will SEE (monoisotopic, not average -- they differ by ~0.66 Da 20// at 1 kDa and that is the difference between a hit and a miss), the m/z of 21// every charge state, and the b/y ladder for identification -- and it 22// REFUSES a sequence containing a residue it does not know rather than 23// silently scoring the unknown as zero mass. 24// 25// PRECISION, stated honestly (not asserted): 26// - Residue masses are rounded to Q4, so an n-mer accumulates up to n/2 27// Q4 units of rounding. Measured against literature: angiotensin II 28// (8-mer) is EXACT (1045.5345); bradykinin (9-mer) lands 1 Q4 unit 29// (0.0001 Da) high. pep_mass_agrees() makes that tolerance explicit. 30// - pI comes from bisection over pep_pow10_q3, whose decade table is 31// 0.1-log granular, so pI is good to about +/-0.05 pH -- adequate to 32// choose a buffer, NOT a substitute for a titration. 33// 34// Grounding (cited; researcher-groundable): 35// iupac_amino_acid_residue_monoisotopic_masses (residue table) 36// emboss_pepstats_pka_set (pKa set for pI/charge) 37// kyte_doolittle_hydropathy_1982 (GRAVY) 38// pace_1995_extinction_coefficient_280nm (Trp 5500 / Tyr 1490 / cystine 125) 39// biemann_1990_peptide_fragment_ion_nomenclature (b/y ions) 40// 41// genealogy_id: peptide_chemistry + nishi_chem_suite + nishi_food_science 42 43import "nx_syscalls.nx" 44import "nx_pow10.nx" 45const PEP_MAGIC_710371: i64 = 710371 46const PEP_MAGIC_1030092: i64 = 1030092 47const PEP_MAGIC_1150269: i64 = 1150269 48const PEP_MAGIC_1290426: i64 = 1290426 49const PEP_MAGIC_1470684: i64 = 1470684 50const PEP_MAGIC_570215: i64 = 570215 51const PEP_MAGIC_1370589: i64 = 1370589 52const PEP_MAGIC_1130841: i64 = 1130841 53const PEP_MAGIC_1280950: i64 = 1280950 54const PEP_MAGIC_1310405: i64 = 1310405 55const PEP_MAGIC_1140429: i64 = 1140429 56const PEP_MAGIC_970528: i64 = 970528 57const PEP_MAGIC_1280586: i64 = 1280586 58const PEP_MAGIC_1561011: i64 = 1561011 59const PEP_MAGIC_870320: i64 = 870320 60const PEP_MAGIC_1010477: i64 = 1010477 61const PEP_MAGIC_990684: i64 = 990684 62const PEP_MAGIC_1860793: i64 = 1860793 63const PEP_MAGIC_1630633: i64 = 1630633 64const PEP_MAGIC_710788: i64 = 710788 65const PEP_MAGIC_1031388: i64 = 1031388 66const PEP_MAGIC_1150886: i64 = 1150886 67const PEP_MAGIC_1291155: i64 = 1291155 68const PEP_MAGIC_1471766: i64 = 1471766 69const PEP_MAGIC_570519: i64 = 570519 70const PEP_MAGIC_1371411: i64 = 1371411 71const PEP_MAGIC_1131594: i64 = 1131594 72const PEP_MAGIC_1281741: i64 = 1281741 73const PEP_MAGIC_1311926: i64 = 1311926 74const PEP_MAGIC_1141038: i64 = 1141038 75const PEP_MAGIC_971167: i64 = 971167 76const PEP_MAGIC_1281307: i64 = 1281307 77const PEP_MAGIC_1561875: i64 = 1561875 78const PEP_MAGIC_870782: i64 = 870782 79const PEP_MAGIC_1011051: i64 = 1011051 80const PEP_MAGIC_991326: i64 = 991326 81const PEP_MAGIC_1862132: i64 = 1862132 82const PEP_MAGIC_1631760: i64 = 1631760 83const PEP_MAGIC_1800: i64 = 1800 84const PEP_MAGIC_2500: i64 = 2500 85const PEP_MAGIC_3500: i64 = 3500 86const PEP_MAGIC_2800: i64 = 2800 87const PEP_MAGIC_3200: i64 = 3200 88const PEP_MAGIC_4500: i64 = 4500 89const PEP_MAGIC_3900: i64 = 3900 90const PEP_MAGIC_3800: i64 = 3800 91const PEP_MAGIC_1900: i64 = 1900 92const PEP_MAGIC_1600: i64 = 1600 93const PEP_MAGIC_4200: i64 = 4200 94const PEP_MAGIC_1300: i64 = 1300 95const PEP_MAGIC_1000000: i64 = 1000000 96 97// ===== Sealed enum: the 20 proteinogenic residues (1-letter order) ===== 98 99const PEP_A: i64 = 0 100const PEP_C: i64 = 1 101const PEP_D: i64 = 2 102const PEP_E: i64 = 3 103const PEP_F: i64 = 4 104const PEP_G: i64 = 5 105const PEP_H: i64 = 6 106const PEP_I: i64 = 7 107const PEP_K: i64 = 8 108const PEP_L: i64 = 9 109const PEP_M: i64 = 10 110const PEP_N: i64 = 11 111const PEP_P: i64 = 12 112const PEP_Q: i64 = 13 113const PEP_R: i64 = 14 114const PEP_S: i64 = 15 115const PEP_T: i64 = 16 116const PEP_V: i64 = 17 117const PEP_W: i64 = 18 118const PEP_Y: i64 = 19 119const PEP_N_AA: i64 = 20 120 121const PEP_INVALID: i64 = 0 - 1 122 123// ===== Physical constants (Q4 micro-AMU) ============================== 124// Water is added once per chain (the residue table is already 125// free-amino-acid MINUS water, so a chain of n residues + 1 water = peptide). 126 127const PEP_WATER_MONO_Q4: i64 = 180106 // H2O monoisotopic 18.010565 Da 128const PEP_WATER_AVG_Q4: i64 = 180153 // H2O average 18.0153 Da 129// Proton mass; same value nx_chem_mass.nx uses for [M+H]+ so ion m/z from 130// a sequence and from a MolGraph are directly comparable. 131const PEP_PROTON_Q4: i64 = 10073 // 1.00728 Da 132 133// ===== Molar extinction at 280 nm (M^-1 cm^-1, Pace 1995) ============= 134 135const PEP_EXT_TRP: i64 = 5500 136const PEP_EXT_TYR: i64 = 1490 137const PEP_EXT_CYSTINE: i64 = 125 // per DISULFIDE, not per free Cys 138 139// ===== pKa set (EMBOSS pepstats), pH x 1000 =========================== 140 141const PEP_PKA_CTERM: i64 = 3600 142const PEP_PKA_NTERM: i64 = 8600 143const PEP_PKA_C: i64 = 8500 144const PEP_PKA_D: i64 = 3900 145const PEP_PKA_E: i64 = 4100 146const PEP_PKA_H: i64 = 6500 147const PEP_PKA_K: i64 = 10800 148const PEP_PKA_R: i64 = 12500 149const PEP_PKA_Y: i64 = 10100 150 151// pI bisection bounds + depth. 20 halvings of a 14-pH span lands far below 152// the 0.1-log granularity of the decade table, so the table is the limit. 153const PEP_PI_LO_MILLI: i64 = 0 154const PEP_PI_HI_MILLI: i64 = 14000 155const PEP_PI_ITERS: i64 = 20 156 157// Beyond +/-6 log units the ionisation fraction is under 1e-6 -- saturate 158// rather than overflow. 159const PEP_POW10_CLAMP: i64 = 6000 160 161// ===== Residue masses ================================================= 162// 163// Monoisotopic residue mass (Q4). Leu and Ile are ISOBARIC (identical 164// mass) -- that is real chemistry, not a copy-paste error: MS cannot 165// distinguish them without fragmentation of the side chain. 166 167func pep_residue_mono_q4(aa: i64) -> i64 { 168 if aa == PEP_A { return PEP_MAGIC_710371 } 169 if aa == PEP_C { return PEP_MAGIC_1030092 } 170 if aa == PEP_D { return PEP_MAGIC_1150269 } 171 if aa == PEP_E { return PEP_MAGIC_1290426 } 172 if aa == PEP_F { return PEP_MAGIC_1470684 } 173 if aa == PEP_G { return PEP_MAGIC_570215 } 174 if aa == PEP_H { return PEP_MAGIC_1370589 } 175 if aa == PEP_I { return PEP_MAGIC_1130841 } 176 if aa == PEP_K { return PEP_MAGIC_1280950 } 177 if aa == PEP_L { return PEP_MAGIC_1130841 } 178 if aa == PEP_M { return PEP_MAGIC_1310405 } 179 if aa == PEP_N { return PEP_MAGIC_1140429 } 180 if aa == PEP_P { return PEP_MAGIC_970528 } 181 if aa == PEP_Q { return PEP_MAGIC_1280586 } 182 if aa == PEP_R { return PEP_MAGIC_1561011 } 183 if aa == PEP_S { return PEP_MAGIC_870320 } 184 if aa == PEP_T { return PEP_MAGIC_1010477 } 185 if aa == PEP_V { return PEP_MAGIC_990684 } 186 if aa == PEP_W { return PEP_MAGIC_1860793 } 187 if aa == PEP_Y { return PEP_MAGIC_1630633 } 188 return 0 189} 190 191// Average (isotope-weighted) residue mass (Q4) -- what a low-resolution 192// instrument or a gravimetric calculation sees. 193func pep_residue_avg_q4(aa: i64) -> i64 { 194 if aa == PEP_A { return PEP_MAGIC_710788 } 195 if aa == PEP_C { return PEP_MAGIC_1031388 } 196 if aa == PEP_D { return PEP_MAGIC_1150886 } 197 if aa == PEP_E { return PEP_MAGIC_1291155 } 198 if aa == PEP_F { return PEP_MAGIC_1471766 } 199 if aa == PEP_G { return PEP_MAGIC_570519 } 200 if aa == PEP_H { return PEP_MAGIC_1371411 } 201 if aa == PEP_I { return PEP_MAGIC_1131594 } 202 if aa == PEP_K { return PEP_MAGIC_1281741 } 203 if aa == PEP_L { return PEP_MAGIC_1131594 } 204 if aa == PEP_M { return PEP_MAGIC_1311926 } 205 if aa == PEP_N { return PEP_MAGIC_1141038 } 206 if aa == PEP_P { return PEP_MAGIC_971167 } 207 if aa == PEP_Q { return PEP_MAGIC_1281307 } 208 if aa == PEP_R { return PEP_MAGIC_1561875 } 209 if aa == PEP_S { return PEP_MAGIC_870782 } 210 if aa == PEP_T { return PEP_MAGIC_1011051 } 211 if aa == PEP_V { return PEP_MAGIC_991326 } 212 if aa == PEP_W { return PEP_MAGIC_1862132 } 213 if aa == PEP_Y { return PEP_MAGIC_1631760 } 214 return 0 215} 216 217// Kyte-Doolittle hydropathy index x 1000 (positive = hydrophobic). 218func pep_hydropathy_q3(aa: i64) -> i64 { 219 if aa == PEP_A { return PEP_MAGIC_1800 } 220 if aa == PEP_C { return PEP_MAGIC_2500 } 221 if aa == PEP_D { return 0 - PEP_MAGIC_3500 } 222 if aa == PEP_E { return 0 - PEP_MAGIC_3500 } 223 if aa == PEP_F { return PEP_MAGIC_2800 } 224 if aa == PEP_G { return 0 - 400 } 225 if aa == PEP_H { return 0 - PEP_MAGIC_3200 } 226 if aa == PEP_I { return PEP_MAGIC_4500 } 227 if aa == PEP_K { return 0 - PEP_MAGIC_3900 } 228 if aa == PEP_L { return PEP_MAGIC_3800 } 229 if aa == PEP_M { return PEP_MAGIC_1900 } 230 if aa == PEP_N { return 0 - PEP_MAGIC_3500 } 231 if aa == PEP_P { return 0 - PEP_MAGIC_1600 } 232 if aa == PEP_Q { return 0 - PEP_MAGIC_3500 } 233 if aa == PEP_R { return 0 - PEP_MAGIC_4500 } 234 if aa == PEP_S { return 0 - 800 } 235 if aa == PEP_T { return 0 - 700 } 236 if aa == PEP_V { return PEP_MAGIC_4200 } 237 if aa == PEP_W { return 0 - 900 } 238 if aa == PEP_Y { return 0 - PEP_MAGIC_1300 } 239 return 0 240} 241 242// ===== Sequence parsing (fail-closed) ================================= 243// 244// An unrecognised letter returns PEP_INVALID rather than mapping to a 245// zero-mass residue -- so a typo or a non-standard residue (B/J/O/U/X/Z, 246// or lowercase) can never fabricate a plausible-looking mass. 247// ASCII codes are used directly (65='A' .. 89='Y'). 248 249func pep_aa_from_char(c: i64) -> i64 { 250 if c == 65 { return PEP_A } 251 if c == 67 { return PEP_C } 252 if c == 68 { return PEP_D } 253 if c == 69 { return PEP_E } 254 if c == 70 { return PEP_F } 255 if c == 71 { return PEP_G } 256 if c == 72 { return PEP_H } 257 if c == 73 { return PEP_I } 258 if c == 75 { return PEP_K } 259 if c == 76 { return PEP_L } 260 if c == 77 { return PEP_M } 261 if c == 78 { return PEP_N } 262 if c == 80 { return PEP_P } 263 if c == 81 { return PEP_Q } 264 if c == 82 { return PEP_R } 265 if c == 83 { return PEP_S } 266 if c == 84 { return PEP_T } 267 if c == 86 { return PEP_V } 268 if c == 87 { return PEP_W } 269 if c == 89 { return PEP_Y } 270 return PEP_INVALID 271} 272 273func pep_len(seq: *u8) -> i64 { 274 var i: i64 = 0 275 while seq[i] != 0 { i = i + 1 } 276 return i 277} 278 279// 1 iff every character maps to a known residue AND the chain is non-empty. 280func pep_seq_valid(seq: *u8) -> i64 { 281 var i: i64 = 0 282 var ok: i64 = 1 283 let n: i64 = pep_len(seq) 284 if n <= 0 { return 0 } 285 while i < n { 286 let c: i64 = seq[i] as i64 287 let aa: i64 = pep_aa_from_char(c) 288 if aa == PEP_INVALID { ok = 0 } 289 i = i + 1 290 } 291 return ok 292} 293 294func pep_count_aa(seq: *u8, aa_want: i64) -> i64 { 295 var i: i64 = 0 296 var k: i64 = 0 297 let n: i64 = pep_len(seq) 298 while i < n { 299 let c: i64 = seq[i] as i64 300 let aa: i64 = pep_aa_from_char(c) 301 if aa == aa_want { k = k + 1 } 302 i = i + 1 303 } 304 return k 305} 306 307// ===== Molecular mass ================================================= 308// 309// Peptide mass = sum(residue masses) + one water (the residues are already 310// water-subtracted, so this is mass conservation over n-1 peptide bonds). 311// REFUSES (returns PEP_INVALID) on an unknown residue. 312 313func pep_mass_mono_q4(seq: *u8) -> i64 { 314 var i: i64 = 0 315 var total: i64 = PEP_WATER_MONO_Q4 316 if pep_seq_valid(seq) != 1 { return PEP_INVALID } 317 let n: i64 = pep_len(seq) 318 while i < n { 319 let c: i64 = seq[i] as i64 320 let aa: i64 = pep_aa_from_char(c) 321 let m: i64 = pep_residue_mono_q4(aa) 322 total = total + m 323 i = i + 1 324 } 325 return total 326} 327 328func pep_mass_avg_q4(seq: *u8) -> i64 { 329 var i: i64 = 0 330 var total: i64 = PEP_WATER_AVG_Q4 331 if pep_seq_valid(seq) != 1 { return PEP_INVALID } 332 let n: i64 = pep_len(seq) 333 while i < n { 334 let c: i64 = seq[i] as i64 335 let aa: i64 = pep_aa_from_char(c) 336 let m: i64 = pep_residue_avg_q4(aa) 337 total = total + m 338 i = i + 1 339 } 340 return total 341} 342 343// Explicit tolerance check against a reference mass -- makes the Q4 344// accumulation budget a TESTABLE claim instead of a promise. 345func pep_mass_agrees(computed_q4: i64, reference_q4: i64, tol_q4: i64) -> i64 { 346 if computed_q4 == PEP_INVALID { return 0 } 347 var d: i64 = computed_q4 - reference_q4 348 if d < 0 { d = 0 - d } 349 if d <= tol_q4 { return 1 } 350 return 0 351} 352 353// ===== Ion m/z ======================================================== 354// 355// [M+nH]n+ : (M + n*proton)/n. Electrospray of a peptide gives a charge 356// LADDER, so the observed m/z depends on n -- this returns the value for 357// any n, which is what deconvoluting a real spectrum needs. 358 359func pep_mz_q4(mass_q4: i64, n_charge: i64) -> i64 { 360 if n_charge <= 0 { return 0 } 361 if mass_q4 == PEP_INVALID { return PEP_INVALID } 362 let protons: i64 = n_charge * PEP_PROTON_Q4 363 let total: i64 = mass_q4 + protons 364 return total / n_charge 365} 366 367// ===== b / y fragment ions (Biemann nomenclature) ===================== 368// 369// Backbone cleavage at the amide bond gives complementary series: 370// b_i = sum of the FIRST i residues + proton (N-terminal) 371// y_i = sum of the LAST i residues + water + proton (C-terminal) 372// Matching an observed peak list against these two ladders IS peptide 373// sequencing; b_i + y_(n-i) = M + 2*proton, which pep_by_complement_ok 374// checks as a self-consistency invariant. 375 376func pep_b_ion_q4(seq: *u8, i_frag: i64) -> i64 { 377 var i: i64 = 0 378 var total: i64 = PEP_PROTON_Q4 379 if pep_seq_valid(seq) != 1 { return PEP_INVALID } 380 let n: i64 = pep_len(seq) 381 if i_frag <= 0 { return PEP_INVALID } 382 if i_frag >= n { return PEP_INVALID } 383 while i < i_frag { 384 let c: i64 = seq[i] as i64 385 let aa: i64 = pep_aa_from_char(c) 386 let m: i64 = pep_residue_mono_q4(aa) 387 total = total + m 388 i = i + 1 389 } 390 return total 391} 392 393func pep_y_ion_q4(seq: *u8, i_frag: i64) -> i64 { 394 var total: i64 = PEP_PROTON_Q4 + PEP_WATER_MONO_Q4 395 if pep_seq_valid(seq) != 1 { return PEP_INVALID } 396 let n: i64 = pep_len(seq) 397 if i_frag <= 0 { return PEP_INVALID } 398 if i_frag >= n { return PEP_INVALID } 399 let start: i64 = n - i_frag 400 var i: i64 = start 401 while i < n { 402 let c: i64 = seq[i] as i64 403 let aa: i64 = pep_aa_from_char(c) 404 let m: i64 = pep_residue_mono_q4(aa) 405 total = total + m 406 i = i + 1 407 } 408 return total 409} 410 411// Invariant: b_i + y_(n-i) == M + 2 protons, for every interior i. 412func pep_by_complement_ok(seq: *u8, i_frag: i64) -> i64 { 413 let n: i64 = pep_len(seq) 414 let b: i64 = pep_b_ion_q4(seq, i_frag) 415 let y: i64 = pep_y_ion_q4(seq, n - i_frag) 416 if b == PEP_INVALID { return 0 } 417 if y == PEP_INVALID { return 0 } 418 let m: i64 = pep_mass_mono_q4(seq) 419 if m == PEP_INVALID { return 0 } 420 let want: i64 = m + 2 * PEP_PROTON_Q4 421 if b + y == want { return 1 } 422 return 0 423} 424 425// ===== Integer decade exponential ===================================== 426// 427// Delegated to nx_pow10 (shared with the thermal-kinetics lane, which needs 428// the same decade-log primitive). The name is kept so callers and the gate 429// are untouched by the extraction. nx_pow10 INTERPOLATES between table 430// entries, so this is now better than the 0.1-log-granular original -- the 431// pI figure is correspondingly tighter than the +/-0.05 pH stated above. 432 433func pep_pow10_q3(d_milli: i64) -> i64 { 434 return ipow10_q3(d_milli) 435} 436 437// ===== Ionisation fractions (Henderson-Hasselbalch) =================== 438// 439// Basic group, fraction PROTONATED (carries +1): 1/(1 + 10^(pH-pKa)) 440// Acidic group, fraction DEPROTONATED (carries -1): 1/(1 + 10^(pKa-pH)) 441// Both returned x1000. At pH == pKa each is exactly 500 (half-ionised), 442// which is the cheapest correctness check there is. 443 444func pep_frac_basic_q3(ph_milli: i64, pka_milli: i64) -> i64 { 445 let p: i64 = pep_pow10_q3(ph_milli - pka_milli) 446 let denom: i64 = 1000 + p 447 if denom <= 0 { return 0 } 448 return PEP_MAGIC_1000000 / denom 449} 450 451func pep_frac_acidic_q3(ph_milli: i64, pka_milli: i64) -> i64 { 452 let p: i64 = pep_pow10_q3(pka_milli - ph_milli) 453 let denom: i64 = 1000 + p 454 if denom <= 0 { return 0 } 455 return PEP_MAGIC_1000000 / denom 456} 457 458// ===== Net charge at a given pH ======================================= 459// 460// Sums the free alpha-amino terminus, the free alpha-carboxyl terminus, 461// and every ionisable side chain. Returned x1000 (so +2500 = +2.5 e). 462 463func pep_charge_q3(seq: *u8, ph_milli: i64) -> i64 { 464 if pep_seq_valid(seq) != 1 { return 0 } 465 let n_k: i64 = pep_count_aa(seq, PEP_K) 466 let n_r: i64 = pep_count_aa(seq, PEP_R) 467 let n_h: i64 = pep_count_aa(seq, PEP_H) 468 let n_d: i64 = pep_count_aa(seq, PEP_D) 469 let n_e: i64 = pep_count_aa(seq, PEP_E) 470 let n_c: i64 = pep_count_aa(seq, PEP_C) 471 let n_y: i64 = pep_count_aa(seq, PEP_Y) 472 var pos: i64 = pep_frac_basic_q3(ph_milli, PEP_PKA_NTERM) 473 let f_k: i64 = pep_frac_basic_q3(ph_milli, PEP_PKA_K) 474 let f_r: i64 = pep_frac_basic_q3(ph_milli, PEP_PKA_R) 475 let f_h: i64 = pep_frac_basic_q3(ph_milli, PEP_PKA_H) 476 pos = pos + n_k * f_k 477 pos = pos + n_r * f_r 478 pos = pos + n_h * f_h 479 var neg: i64 = pep_frac_acidic_q3(ph_milli, PEP_PKA_CTERM) 480 let f_d: i64 = pep_frac_acidic_q3(ph_milli, PEP_PKA_D) 481 let f_e: i64 = pep_frac_acidic_q3(ph_milli, PEP_PKA_E) 482 let f_c: i64 = pep_frac_acidic_q3(ph_milli, PEP_PKA_C) 483 let f_y: i64 = pep_frac_acidic_q3(ph_milli, PEP_PKA_Y) 484 neg = neg + n_d * f_d 485 neg = neg + n_e * f_e 486 neg = neg + n_c * f_c 487 neg = neg + n_y * f_y 488 return pos - neg 489} 490 491// ===== Isoelectric point ============================================== 492// 493// pI = the pH where net charge crosses zero. Charge is monotonically 494// DECREASING in pH (every group deprotonates as pH rises), so plain 495// bisection is sound -- no derivative, no float, no iteration limit risk. 496 497func pep_pi_milli(seq: *u8) -> i64 { 498 if pep_seq_valid(seq) != 1 { return PEP_INVALID } 499 var lo: i64 = PEP_PI_LO_MILLI 500 var hi: i64 = PEP_PI_HI_MILLI 501 var i: i64 = 0 502 while i < PEP_PI_ITERS { 503 let mid: i64 = (lo + hi) / 2 504 let c: i64 = pep_charge_q3(seq, mid) 505 if c > 0 { lo = mid } 506 if c <= 0 { hi = mid } 507 i = i + 1 508 } 509 return (lo + hi) / 2 510} 511 512// ===== UV absorbance at 280 nm ======================================== 513// 514// Molar extinction is EXACT integer arithmetic (Pace 1995): only Trp, Tyr 515// and disulfide-bonded cystine absorb. n_cystine is a DISULFIDE count and 516// must be supplied by the caller -- it is not derivable from sequence 517// alone, and guessing it would be fabrication. A peptide with no W and no 518// Y has extinction 0, meaning A280 CANNOT be used to quantify it: a real 519// and frequently-missed laboratory constraint. 520 521func pep_ext_coeff_280(seq: *u8, n_cystine: i64) -> i64 { 522 if pep_seq_valid(seq) != 1 { return PEP_INVALID } 523 let n_w: i64 = pep_count_aa(seq, PEP_W) 524 let n_y: i64 = pep_count_aa(seq, PEP_Y) 525 var e: i64 = n_w * PEP_EXT_TRP 526 e = e + n_y * PEP_EXT_TYR 527 if n_cystine > 0 { e = e + n_cystine * PEP_EXT_CYSTINE } 528 return e 529} 530 531// 1 iff A280 is a usable quantitation method for this peptide. 532func pep_a280_quantifiable(seq: *u8, n_cystine: i64) -> i64 { 533 let e: i64 = pep_ext_coeff_280(seq, n_cystine) 534 if e == PEP_INVALID { return 0 } 535 if e > 0 { return 1 } 536 return 0 537} 538 539// ===== GRAVY (grand average of hydropathy) ============================ 540// 541// Mean Kyte-Doolittle index over the chain, x1000. Positive = hydrophobic 542// (aggregation / membrane-association risk); negative = soluble. 543 544func pep_gravy_q3(seq: *u8) -> i64 { 545 var i: i64 = 0 546 var sum: i64 = 0 547 if pep_seq_valid(seq) != 1 { return PEP_INVALID } 548 let n: i64 = pep_len(seq) 549 if n <= 0 { return PEP_INVALID } 550 while i < n { 551 let c: i64 = seq[i] as i64 552 let aa: i64 = pep_aa_from_char(c) 553 let h: i64 = pep_hydropathy_q3(aa) 554 sum = sum + h 555 i = i + 1 556 } 557 return sum / n 558} 559 560func pep_is_hydrophobic(seq: *u8) -> i64 { 561 let g: i64 = pep_gravy_q3(seq) 562 if g == PEP_INVALID { return 0 } 563 if g > 0 { return 1 } 564 return 0 565}