OSCR

Cryo-EM insights into isoform-specific properties of the IP<sub>3</sub>R2 channel.

Code ↔ Paper

2 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 2 matches
  1. [1] § Methods › Cryo-EM data acquisition, image processing, and model building ↔ pyHole.py, lines 137–145 · score 0.53 · van der Waals, radius
  2. [2] § Methods › Cryo-EM data acquisition, image processing, and model building ↔ pyhole_chimerax_fix4/src/holepy.py, lines 139–147 · score 0.53 · van der Waals, radius

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Python · 1,089 lines · 48 KB · GPL-3.0 · 1 match

  1. """pyHole: HOLE-like pore profile calculation with straight/curved centerline options.
  3. Standalone command-line tool. Accepts legacy PDB directly, or any other
  4. gemmi-readable structure format (mmCIF, mmJSON, ...), which is transparently
  5. converted to a temp legacy PDB file before parsing.
  6. """
  7. from pathlib import Path
  8. import argparse
  9. import sys
  10. import re
  11. import math
  12. import json
  13. import csv
  14. import os
  15. import tempfile
  16. from typing import Dict, List, Optional, Tuple
  17. import numpy as np
  18. from scipy.spatial import cKDTree
  19. from scipy.ndimage import label as _cc_label
  20. import gemmi
  21. import mrcfile
  23. __version__ = "1.1"
  25. # --- Constants (in sync with CLI v8) ---
  26. ONE_LETTER = {
  27. 'ALA': 'A', 'CYS': 'C', 'ASP': 'D', 'GLU': 'E', 'PHE': 'F', 'GLY': 'G',
  28. 'HIS': 'H', 'ILE': 'I', 'LYS': 'K', 'LEU': 'L', 'MET': 'M', 'ASN': 'N',
  29. 'PRO': 'P', 'GLN': 'Q', 'ARG': 'R', 'SER': 'S', 'THR': 'T', 'VAL': 'V',
  30. 'TRP': 'W', 'TYR': 'Y', 'SEC': 'U', 'PYL': 'O',
  31. }
  32. VDW_DEFAULT = {
  33. 'H': 1.20, 'C': 1.70, 'N': 1.55, 'O': 1.52, 'F': 1.47, 'P': 1.80,
  34. 'S': 1.80, 'CL': 1.75, 'BR': 1.85, 'I': 1.98, 'NA': 2.27, 'MG': 1.73,
  35. 'K': 2.75, 'CA': 2.31, 'ZN': 1.39, 'FE': 1.56, 'CU': 1.40, 'MN': 1.61,
  36. }
  37. KD_HYDRO = {
  38. 'ILE': 4.5, 'VAL': 4.2, 'LEU': 3.8, 'PHE': 2.8, 'CYS': 2.5, 'MET': 1.9,
  39. 'ALA': 1.8, 'GLY': -0.4, 'THR': -0.7, 'SER': -0.8, 'TRP': -0.9,
  40. 'TYR': -1.3, 'PRO': -1.6, 'HIS': -3.2, 'GLU': -3.5, 'GLN': -3.5,
  41. 'ASP': -3.5, 'ASN': -3.5, 'LYS': -3.9, 'ARG': -4.5,
  42. }
  43. KD_MIN = min(KD_HYDRO.values())
  44. KD_MAX = max(KD_HYDRO.values())
  45. CHARGE_MAP = {'ARG': +1.0, 'LYS': +1.0, 'HIS': +0.1, 'ASP': -1.0, 'GLU': -1.0}
  46. ELEC_MIN, ELEC_MAX = -1.0, +1.0 # mean formal charge per slice maps into [-1, 1] before optional 01 remap
  47. # Floor for reported mesh/centerline radius (Å); avoids zero/negative clearance breaking mesh / coloring.
  48. MIN_PORE_RADIUS_A = 0.1
  51. # --- mmCIF / any-gemmi-format -> legacy PDB (self-contained; see
  52. # cryomodel/pore/structure_io.py for the full rationale) ---
  53. _MAX_ATOMS = 99_999
  54. _MAX_RESI = 9_999
  55. _MIN_RESI = -999
  56. _LEGACY_PDB_SUFFIXES = {'.pdb', '.ent', '.pdb1'}
  59. class StructureTooLargeForLegacyPDB(ValueError):
  60. """Raised when a structure can't be safely represented as legacy PDB."""
  63. def ensure_legacy_pdb(path, tmp_dir: Optional[str] = None) -> Path:
  64. """Return a path to a legacy fixed-width PDB file equivalent to ``path``.
  66. holepy.py's parser assumes fixed-width legacy-PDB columns and does not
  67. understand hybrid-36 encoding, so structures that would overflow legacy
  68. PDB's fields (>99,999 atoms, >4-digit residue numbers, multi-character
  69. chain IDs) raise StructureTooLargeForLegacyPDB instead of silently
  70. corrupting atom/residue numbering. Subset the structure and retry.
  71. """
  72. path = Path(path)
  73. if path.suffix.lower() in _LEGACY_PDB_SUFFIXES:
  74. return path
  76. st = gemmi.read_structure(str(path))
  77. if len(st) == 0:
  78. raise ValueError(f"No models found in '{path.name}'.")
  80. single = gemmi.Structure()
  81. single.name = st.name
  82. single.cell = st.cell
  83. single.spacegroup_hm = st.spacegroup_hm
  84. single.add_model(st[0])
  85. single.setup_entities()
  87. n_atoms = 0
  88. bad_chains = set()
  89. bad_residues = []
  90. for chain in single[0]:
  91. if len(chain.name) > 1:
  92. bad_chains.add(chain.name)
  93. for res in chain:
  94. if res.seqid.num > _MAX_RESI or res.seqid.num < _MIN_RESI:
  95. bad_residues.append((chain.name, res.seqid.num))
  96. n_atoms += len(res)
  98. problems = []
  99. if n_atoms > _MAX_ATOMS:
  100. problems.append(f"{n_atoms} atoms (legacy PDB's serial field caps at {_MAX_ATOMS})")
  101. if bad_chains:
  102. problems.append(f"chain ID(s) longer than 1 character: {sorted(bad_chains)}")
  103. if bad_residues:
  104. sample = ', '.join(f"{c}:{n}" for c, n in bad_residues[:5])
  105. more = f" (+{len(bad_residues) - 5} more)" if len(bad_residues) > 5 else ""
  106. problems.append(f"{len(bad_residues)} residue number(s) outside PDB's 4-digit field, e.g. {sample}{more}")
  108. if problems:
  109. raise StructureTooLargeForLegacyPDB(
  110. f"'{path.name}' cannot be safely converted to legacy PDB for pyHole: " + "; ".join(problems) + ". "
  111. "Subset the structure (e.g. to the chain(s)/region around your pore) and retry."
  112. )
  114. fd, tmp_path = tempfile.mkstemp(suffix='.pdb', prefix='pyhole_cif2pdb_', dir=tmp_dir)
  115. os.close(fd)
  116. single.write_pdb(tmp_path)
  117. return Path(tmp_path)
  120. def guess_element(atom_name, element_field):
  121. if element_field and element_field.strip():
  122. return element_field.strip().upper()
  123. an = atom_name.strip()
  124. if not an:
  125. return 'C'
  126. if an[0].isdigit() and len(an) >= 2:
  127. c = an[1]
  128. if len(an) >= 3 and an[2].isalpha():
  129. return (c + an[2]).upper()
  130. return c.upper()
  131. c = an[0]
  132. if len(an) >= 2 and an[1].isalpha() and an[1].islower():
  133. return (c + an[1]).upper()
  134. return c.upper()
  137. def vdw_radius(element, custom):
  138. """Get van der Waals radius for element.
  140. ``custom`` must already have upper-cased keys -- callers should normalize
  141. a loaded JSON file with normalize_vdw_keys() before passing it in here,
  142. since this does a case-sensitive lookup after upper-casing ``element``.
  143. """
  144. e = element.upper()
  145. return custom.get(e, VDW_DEFAULT.get(e, 1.7))
  148. def normalize_vdw_keys(raw: Dict[str, float]) -> Dict[str, float]:
  149. """Upper-case all keys in a custom VDW radii dict.
  151. Without this, a JSON file with lowercase element keys (e.g. {"na": 1.1})
  152. silently fails every lookup in vdw_radius() and falls back to the
  153. generic default -- the custom radii are ignored with no warning.
  154. """
  155. return {str(k).upper(): v for k, v in (raw or {}).items()}
  158. class Atom:
  159. __slots__ = ('x', 'y', 'z', 'name', 'resname', 'chain', 'resi', 'icode', 'element', 'occ')
  161. def __init__(self, x, y, z, name, resname, chain, resi, icode, element, occ):
  162. self.x = float(x)
  163. self.y = float(y)
  164. self.z = float(z)
  165. self.name = name.strip()
  166. self.resname = resname
  167. self.chain = chain
  168. self.resi = int(resi)
  169. self.icode = icode # PDB insertion code (single char, '' if none)
  170. self.element = element
  171. self.occ = float(occ)
  174. def load_pdb_atoms(path, include_h=True, include_hetatm=True):
  175. atoms = []
  176. with open(path, 'r') as f:
  177. for line in f:
  178. if line[:6] not in ('ATOM ', 'HETATM'):
  179. continue
  180. if (not include_hetatm) and line[:6] == 'HETATM':
  181. continue
  182. name = line[12:16]
  183. resname = line[17:20].strip()
  184. chain = line[21].strip() or 'A'
  185. resi_str = line[22:26].strip() or '0'
  186. # Column 27 (index 26) is the PDB insertion code. Residues that
  187. # differ only by insertion code (e.g. 55 and 55A) are distinct
  188. # residues -- dropping this silently merged them.
  189. icode = line[26].strip() if len(line) > 26 else ''
  190. try:
  191. resi = int(resi_str)
  192. except ValueError:
  193. continue
  194. x = float(line[30:38])
  195. y = float(line[38:46])
  196. z = float(line[46:54])
  197. occ = float(line[54:60]) if line[54:60].strip() else 1.0
  198. elem_field = line[76:78] if len(line) >= 78 else ''
  199. elem = guess_element(name, elem_field)
  200. if (not include_h) and elem.upper() == 'H':
  201. continue
  202. atoms.append(Atom(x, y, z, name, resname, chain, resi, icode, elem, occ))
  203. return atoms
  206. def parse_residue_tokens(s):
  207. tokens = []
  208. if not s:
  209. return tokens
  210. for t in s.split(','):
  211. t = t.strip()
  212. if not t:
  213. continue
  214. # ONELETTER...NUM/CHAIN or THREELETNUM/CHAIN (chain optional -> match any chain_id)
  215. m = re.match(r'^([A-Za-z]{1,3})(\d+)\s*/\s*([A-Za-z0-9]*)$', t)
  216. if m:
  217. chain = m.group(3).strip() or '*'
  218. tokens.append((chain, int(m.group(2))))
  219. continue
  220. m = re.match(r'^([A-Za-z0-9])\s*[:\s]\s*(\d+)$', t)
  221. if m:
  222. tokens.append((m.group(1), int(m.group(2))))
  223. continue
  224. m = re.match(r'^\s*(\d+)\s*$', t)
  225. if m:
  226. tokens.append(('*', int(m.group(1))))
  227. continue
  228. raise ValueError(f"Could not parse residue token: '{t}'.")
  229. return tokens
  232. def ca_positions_for(atoms, sel):
  233. out = []
  234. for (ch, rnum) in sel:
  235. for a in atoms:
  236. if a.name == 'CA' and a.resi == rnum and (ch == '*' or a.chain == ch):
  237. out.append([a.x, a.y, a.z])
  238. return np.array(out, float) if out else np.zeros((0, 3), float)
  241. def orthonormal_basis_from_axis(axis):
  242. u = axis / (np.linalg.norm(axis) + 1e-12)
  243. a = np.array([1.0, 0.0, 0.0]) if abs(u[0]) < 0.9 else np.array([0.0, 1.0, 0.0])
  244. v = np.cross(u, a)
  245. n = np.linalg.norm(v)
  246. if n < 1e-8:
  247. a = np.array([0.0, 0.0, 1.0])
  248. v = np.cross(u, a)
  249. n = np.linalg.norm(v)
  250. v /= (n + 1e-12)
  251. w = np.cross(u, v)
  252. w /= (np.linalg.norm(w) + 1e-12)
  253. return v, w
  256. class SpatialIndex:
  257. """KD-tree wrapper for fast local clearance queries.
  259. Evaluating a slice used to scan every atom in the structure (O(N) per
  260. query point), which dominates runtime for large assemblies with many
  261. query points (dense/adaptive sampling, curved-centerline refinement).
  262. We instead query only atoms within a radius of the candidate point and
  263. grow that radius until we can prove no atom outside the search ball
  264. could produce a smaller clearance -- returns exactly the same value as
  265. the brute-force scan, just faster.
  266. """
  268. def __init__(self, atom_xyz: np.ndarray, atom_r: np.ndarray, initial_radius: float = 12.0):
  269. self.atom_xyz = atom_xyz
  270. self.atom_r = atom_r
  271. self.tree = cKDTree(atom_xyz) if len(atom_xyz) else None
  272. self.initial_radius = initial_radius
  273. self.max_atom_r = float(np.max(atom_r)) if len(atom_r) else 0.0
  275. def clearance_and_mask(self, c: np.ndarray, contact_eps: float = 0.0):
  276. if self.tree is None or len(self.atom_xyz) == 0:
  277. return float('inf'), np.zeros(0, dtype=bool)
  278. r = self.initial_radius
  279. n_atoms = len(self.atom_xyz)
  280. while True:
  281. idx = self.tree.query_ball_point(c, r=r)
  282. if idx:
  283. idx = np.asarray(idx)
  284. dv = self.atom_xyz[idx] - c
  285. d = np.sqrt((dv * dv).sum(axis=1)) - self.atom_r[idx]
  286. d_min = float(np.min(d))
  287. # Any atom outside this ball is at distance > r, so its
  288. # clearance is > r - max_atom_r. If our current best is
  289. # already <= that bound, no unseen atom can beat it.
  290. if d_min <= (r - self.max_atom_r) or len(idx) >= n_atoms:
  291. full_d = np.full(n_atoms, np.inf)
  292. full_d[idx] = d
  293. mask = full_d <= (d_min + contact_eps)
  294. return d_min, mask
  295. r *= 1.5
  297. def batch_knn_clearance(self, points: np.ndarray, k: int = 12) -> np.ndarray:
  298. """Vectorized clearance (dist - vdw_r) to the nearest atom for many
  299. points at once, via a k-nearest-neighbor search rather than the exact
  300. expanding-radius search used by clearance_and_mask().
  302. This is an approximation: it's only exact if the atom that actually
  303. minimizes clearance is among the k nearest by *raw* distance, which
  304. holds in practice because VDW radii vary within a bounded range
  305. (~1-2.8 A) -- an atom much farther away by raw distance essentially
  306. never wins on clearance. Used for area/volume rasterization and the
  307. density-map surface, where we evaluate thousands-to-millions of grid
  308. points and the expanding-radius search's per-point Python loop would
  309. be too slow; not used for the single-point clearance that determines
  310. the reported minimum pore radius (clearance_and_mask stays exact there).
  311. """
  312. if self.tree is None or len(self.atom_xyz) == 0:
  313. return np.full(len(points), np.inf)
  314. k_eff = min(k, len(self.atom_r))
  315. dists, idxs = self.tree.query(points, k=k_eff)
  316. if k_eff == 1:
  317. dists = dists[:, None]
  318. idxs = idxs[:, None]
  319. clearances = dists - self.atom_r[idxs]
  320. return clearances.min(axis=1)
  323. def _clearance_at_point(atom_xyz, atom_r, c, index: Optional[SpatialIndex] = None):
  324. if index is not None:
  325. d_min, _ = index.clearance_and_mask(c)
  326. return d_min
  327. dv = atom_xyz - c
  328. d = np.sqrt((dv * dv).sum(axis=1)) - atom_r
  329. return float(np.min(d))
  332. def _evaluate_slice(atom_xyz, atom_r, atom_meta, c, contact_eps, hydro_scale, electro_scale,
  333. index: Optional[SpatialIndex] = None):
  334. if index is not None:
  335. r_raw, mask = index.clearance_and_mask(c, contact_eps)
  336. else:
  337. dv = atom_xyz - c
  338. d = np.sqrt((dv * dv).sum(axis=1)) - atom_r
  339. r_raw = float(np.min(d))
  340. mask = d <= (r_raw + contact_eps)
  342. seen = set()
  343. tags = []
  344. hyd = []
  345. chg = []
  346. for i, ok in enumerate(mask):
  347. if not ok:
  348. continue
  349. chain, resname, resi, icode = atom_meta[i]
  350. # Insertion code is part of residue identity: 55 and 55A must not
  351. # be collapsed into a single contributor tag.
  352. key = (chain, resi, icode, resname)
  353. if key in seen:
  354. continue
  355. seen.add(key)
  356. resi_label = f"{resi}{icode}" if icode else str(resi)
  357. tags.append(f"{ONE_LETTER.get(resname.upper(), '?')}{resi_label}/{chain}")
  358. hyd.append(KD_HYDRO.get(resname.upper(), 0.0))
  359. chg.append(CHARGE_MAP.get(resname.upper(), 0.0))
  361. hydro = float(np.mean(hyd)) if hyd else 0.0
  362. electro = float(np.mean(chg)) if chg else 0.0
  363. if hydro_scale == '01':
  364. hydro = (hydro - KD_MIN) / (KD_MAX - KD_MIN) if KD_MAX > KD_MIN else 0.0
  365. if electro_scale == '01':
  366. electro = (electro - ELEC_MIN) / (ELEC_MAX - ELEC_MIN) if ELEC_MAX > ELEC_MIN else 0.5
  368. rmin = max(r_raw, MIN_PORE_RADIUS_A)
  369. return rmin, ';'.join(tags[:50]), hydro, electro
  372. def _trapz(y, x) -> float:
  373. """Trapezoidal integral of y(x). Hand-rolled instead of np.trapz/np.trapezoid
  374. since that function's name changed between numpy versions."""
  375. y = np.asarray(y, dtype=float)
  376. x = np.asarray(x, dtype=float)
  377. if len(y) < 2:
  378. return 0.0
  379. return float(np.sum((y[1:] + y[:-1]) * 0.5 * (x[1:] - x[:-1])))
  382. def slice_area(index: 'SpatialIndex', center: np.ndarray, v: np.ndarray, w: np.ndarray,
  383. probe: float = 0.0, seed_clearance: Optional[float] = None,
  384. min_half_extent: float = 4.0, max_half_extent: float = 20.0,
  385. radius_factor: float = 3.0, grid_step: float = 0.25, k: int = 12):
  386. """Rasterized true cross-sectional area of the pore in the plane through
  387. ``center`` spanned by orthonormal (v, w), instead of treating the slice
  388. as a circle of the single reported min-clearance radius.
  390. A single "radius" number assumes a circular cross-section; real pores are
  391. often lobed, eccentric, or locally widen off-axis, which a circular
  392. assumption hides entirely and which HOLE-style single-radius profiles
  393. can't represent. This grids the plane, marks each cell as accessible
  394. (clearance > probe) or not, and keeps only the connected component
  395. reachable from ``center`` by flood fill -- so an unrelated cavity or
  396. pocket that happens to intersect the same cutting plane elsewhere in the
  397. structure isn't counted as part of the pore.
  399. Returns (area_geometric_A2, area_accessible_A2, half_extent_used_A).
  400. """
  401. if seed_clearance is None:
  402. seed_clearance = float(index.batch_knn_clearance(center[None, :], k=k)[0])
  403. half_extent = float(np.clip(radius_factor * max(seed_clearance, 0.5), min_half_extent, max_half_extent))
  405. n = max(3, int(round(2 * half_extent / grid_step)) + 1)
  406. axis_vals = np.linspace(-half_extent, half_extent, n)
  407. gv, gw = np.meshgrid(axis_vals, axis_vals, indexing='ij')
  408. pts = center[None, None, :] + gv[..., None] * v[None, None, :] + gw[..., None] * w[None, None, :]
  409. clearance = index.batch_knn_clearance(pts.reshape(-1, 3), k=k).reshape(n, n)
  411. cell_area = grid_step * grid_step
  412. center_idx = (n // 2, n // 2) # axis_vals is symmetric about 0 -> this cell is (0, 0)
  414. def flood_area(mask):
  415. if not mask[center_idx]:
  416. # Center point comes from the already-validated centerline, so
  417. # this shouldn't normally trigger; fail safe to 0 rather than
  418. # grabbing an unrelated component if it ever does.
  419. return 0.0
  420. labeled, _ = _cc_label(mask, structure=np.ones((3, 3)))
  421. comp_id = labeled[center_idx]
  422. return float(np.sum(labeled == comp_id)) * cell_area if comp_id else 0.0
  424. area_geom = flood_area(clearance > 0.0)
  425. area_access = flood_area(clearance > probe) if probe > 0 else area_geom
  426. return area_geom, area_access, half_extent
  429. def build_density_grid(index: 'SpatialIndex', rows: List[Dict],
  430. probe: float = 0.0, voxel_size: float = 0.5, margin: float = 3.0,
  431. min_half_extent: float = 4.0, max_half_extent: float = 20.0,
  432. radius_factor: float = 1.3, max_voxels: int = 40_000_000, k: int = 12):
  433. """Build a 3D scalar field (clearance-to-nearest-atom minus probe radius)
  434. over a padded box around the sampled centerline, for isosurfacing the
  435. true pore lumen (level=0 is exactly the probe-accessible boundary)
  436. instead of the tube-of-spheres mesh approximation.
  438. Masked to a swept tube around the centerline: any grid point farther from
  439. its nearest centerline sample than that sample's local half-extent
  440. (+ margin) is set to a large negative sentinel value instead of its raw
  441. clearance. Without this, isosurfacing raw clearance over the *entire*
  442. padded box picks up every ordinary side-chain/helix packing gap in the
  443. surrounding protein, not just the pore -- producing a "swiss cheese"
  444. density full of unrelated cavities rather than one clean pore-shaped
  445. tube. `radius_factor` here defaults smaller/more conservative than the
  446. one used for the 2D area rasterization (slice_area's area_radius_factor,
  447. default 3.0): a too-generous tube here re-admits the same packing-gap
  448. noise this masking exists to remove.
  450. Returns (values[nx,ny,nz] as (z,y,x)-ordered for MRC, origin_xyz, voxel_size).
  451. """
  452. centers = np.array([[r['x'], r['y'], r['z']] for r in rows], dtype=float)
  453. # Per-slice half-extent estimate (same adaptive sizing as slice_area) sets
  454. # both the swept-tube radius and how far sideways from the centerline the
  455. # grid needs to cover.
  456. seed_clear = index.batch_knn_clearance(centers, k=k)
  457. half_extents = np.clip(radius_factor * np.maximum(seed_clear, 0.5), min_half_extent, max_half_extent)
  458. pad = float(np.max(half_extents)) + margin
  460. lo = centers.min(axis=0) - pad
  461. hi = centers.max(axis=0) + pad
  463. nx = max(2, int(np.ceil((hi[0] - lo[0]) / voxel_size)) + 1)
  464. ny = max(2, int(np.ceil((hi[1] - lo[1]) / voxel_size)) + 1)
  465. nz = max(2, int(np.ceil((hi[2] - lo[2]) / voxel_size)) + 1)
  466. n_total = nx * ny * nz
  467. if n_total > max_voxels:
  468. raise ValueError(
  469. f"Requested density grid is {nx}x{ny}x{nz} = {n_total:,} voxels, over the "
  470. f"{max_voxels:,} safety limit. Increase --surface-voxel (coarser spacing) or "
  471. f"reduce --surface-margin / --surface-radius-factor."
  472. )
  474. xs = lo[0] + voxel_size * np.arange(nx)
  475. ys = lo[1] + voxel_size * np.arange(ny)
  476. zs = lo[2] + voxel_size * np.arange(nz)
  477. gx, gy, gz = np.meshgrid(xs, ys, zs, indexing='ij')
  478. pts = np.stack([gx.ravel(), gy.ravel(), gz.ravel()], axis=1)
  480. clearance = index.batch_knn_clearance(pts, k=k)
  482. # Swept-tube mask: find each grid point's nearest centerline sample and
  483. # that sample's local half-extent, and blank out anything farther away
  484. # than that (+ margin) so isosurfacing can't pick up unrelated cavities
  485. # elsewhere in the structure that happen to fall inside the padded box.
  486. centerline_tree = cKDTree(centers)
  487. dist_to_centerline, nearest_idx = centerline_tree.query(pts, k=1)
  488. allowed_radius = half_extents[nearest_idx] + margin
  489. outside_tube = dist_to_centerline > allowed_radius
  491. SENTINEL = -1000.0
  492. values = (clearance - probe).astype(np.float32)
  493. values[outside_tube] = SENTINEL
  494. values = values.reshape(nx, ny, nz)
  495. values_zyx = np.transpose(values, (2, 1, 0)) # mrcfile wants (z, y, x)
  496. return values_zyx, lo, voxel_size
  499. def write_density_map(path, origin: np.ndarray, voxel_size: float, values_zyx: np.ndarray):
  500. """Write a scalar field to a non-crystallographic MRC map (MRC2014 ORIGIN
  501. convention), so it opens already aligned with the source structure in
  502. ChimeraX/ChimeraX-compatible viewers. level=0 is the probe-accessible
  503. pore boundary by construction (see build_density_grid)."""
  504. with mrcfile.new(str(path), overwrite=True) as mrc:
  505. mrc.set_data(np.ascontiguousarray(values_zyx, dtype=np.float32))
  506. mrc.voxel_size = float(voxel_size)
  507. mrc.header.origin.x = float(origin[0])
  508. mrc.header.origin.y = float(origin[1])
  509. mrc.header.origin.z = float(origin[2])
  510. mrc.update_header_stats()
  513. def profile_along_axis(atom_xyz, atom_r, c0, c1, step, eps, atom_meta,
  514. adaptive=False, slope_thresh=0.5, max_refine=3,
  515. hydro_scale='raw', electro_scale='raw', occupancy_metric='hydro',
  516. true_area=True, probe=0.0, area_grid_step=0.25,
  517. area_min_radius=4.0, area_max_radius=20.0, area_radius_factor=3.0):
  518. axis = c1 - c0
  519. L = np.linalg.norm(axis)
  520. if L < 1e-6:
  521. raise ValueError("Top and bottom centers are too close.")
  522. u = axis / L
  523. index = SpatialIndex(atom_xyz, atom_r)
  524. v_global, w_global = orthonormal_basis_from_axis(u)
  526. def ctr(s):
  527. return c0 + u * s
  529. svals = list(np.linspace(0.0, L, max(1, int(round(L / step)) + 1)))
  531. def eval_rows(vals):
  532. out = []
  533. for s in vals:
  534. c = ctr(s)
  535. rmin, tags, hyd, elec = _evaluate_slice(
  536. atom_xyz, atom_r, atom_meta, c, eps, hydro_scale, electro_scale, index=index)
  537. out.append({
  538. 's_A': float(s), 'x': float(c[0]), 'y': float(c[1]), 'z': float(c[2]),
  539. 'radius_A': float(rmin), 'hydro_index': float(hyd), 'electro_index': float(elec),
  540. 'contributors': tags,
  541. })
  542. return out
  544. rows = eval_rows(svals)
  545. if adaptive:
  546. for _ in range(max_refine):
  547. rows.sort(key=lambda r: r['s_A'])
  548. ns = []
  549. for i in range(len(rows) - 1):
  550. s0, r0 = rows[i]['s_A'], rows[i]['radius_A']
  551. s1, r1 = rows[i + 1]['s_A'], rows[i + 1]['radius_A']
  552. ds = s1 - s0
  553. if ds <= step / 2:
  554. continue
  555. slope = abs((r1 - r0) / ds) if ds > 1e-9 else 0.0
  556. if slope > slope_thresh:
  557. ns.append(0.5 * (s0 + s1))
  558. if not ns:
  559. break
  560. rows += eval_rows(ns)
  561. rows.sort(key=lambda r: r['s_A'])
  562. if true_area:
  563. for r in rows:
  564. c = np.array([r['x'], r['y'], r['z']], float)
  565. area_geom, area_access, half_extent = slice_area(
  566. index, c, v_global, w_global, probe=probe,
  567. min_half_extent=area_min_radius, max_half_extent=area_max_radius,
  568. radius_factor=area_radius_factor, grid_step=area_grid_step)
  569. r['area_A2'] = float(area_geom)
  570. r['area_access_A2'] = float(area_access)
  571. r['area_half_extent_A'] = float(half_extent)
  572. for r in rows:
  573. if occupancy_metric == 'hydro':
  574. r['occ_value'] = float(r['hydro_index'])
  575. elif occupancy_metric == 'electro':
  576. r['occ_value'] = float(r['electro_index'])
  577. else:
  578. r['occ_value'] = float(r['radius_A'])
  579. return rows, u, L
  582. def construct_centers_curved(atom_xyz, atom_r, c0, c1, step, curve_radius, curve_iters):
  583. axis = c1 - c0
  584. L = float(np.linalg.norm(axis))
  585. if L < 1e-6:
  586. raise ValueError("Top and bottom centers are too close.")
  587. u = axis / L
  588. v, w = orthonormal_basis_from_axis(u)
  589. index = SpatialIndex(atom_xyz, atom_r)
  590. svals = list(np.linspace(0.0, L, max(1, int(round(L / step)) + 1)))
  591. centers = []
  592. c_prev = c0.copy()
  593. for idx, s in enumerate(svals):
  594. if s <= 1e-9:
  595. c = c0.copy()
  596. elif abs(s - L) <= 1e-9:
  597. c = c1.copy()
  598. else:
  599. ds = svals[idx] - svals[idx - 1]
  600. c = c_prev + u * ds
  601. r = float(curve_radius)
  602. for _ in range(int(curve_iters)):
  603. best_c = c
  604. best_cl = _clearance_at_point(atom_xyz, atom_r, c, index=index)
  605. for dx, dy in [(1, 0), (-1, 0), (0, 1), (0, -1), (1, 1), (-1, 1), (1, -1), (-1, -1), (0, 0)]:
  606. cand = c + (dx * r) * v + (dy * r) * w
  607. cl = _clearance_at_point(atom_xyz, atom_r, cand, index=index)
  608. if cl > best_cl:
  609. best_cl = cl
  610. best_c = cand
  611. c = best_c
  612. r *= 0.5
  613. centers.append(c)
  614. c_prev = c
  615. return centers, u, L
  618. def sample_straight_centers(c0: np.ndarray, c1: np.ndarray, step: float) -> List[np.ndarray]:
  619. """Evenly-spaced points along the straight line from c0 to c1 (inclusive
  620. of both ends). This is the same axis-sampling profile_along_axis does
  621. internally, pulled out standalone so a straight segment can be one leg
  622. of a multi-waypoint path (see build_waypoint_centers)."""
  623. axis = c1 - c0
  624. L = float(np.linalg.norm(axis))
  625. if L < 1e-6:
  626. raise ValueError("Segment endpoints are too close together.")
  627. u = axis / L
  628. svals = np.linspace(0.0, L, max(1, int(round(L / step)) + 1))
  629. return [c0 + u * s for s in svals]
  632. def build_waypoint_centers(atom_xyz, atom_r, waypoints: List[np.ndarray], step: float,
  633. centerline: str, curve_radius: float, curve_iters: float) -> List[np.ndarray]:
  634. """Join 2+ waypoints (bottom, [mid ...], top) into one continuous list of
  635. centerline points, sampling each consecutive pair as either a straight
  636. segment or a curved (locally-refined) segment depending on `centerline`.
  638. This is the "midpoint" feature: an optional intermediate waypoint lets a
  639. user who already knows roughly where a channel bends route the profile
  640. through it, without pyHole having to *find* the bend itself the way
  641. --centerline explore's search does. Joining is genuinely just
  642. concatenation -- arc length (s_A) and local tangent direction are
  643. already computed generically over an arbitrary point sequence by
  644. profile_along_centers, so nothing downstream needs to know a path was
  645. assembled from more than one segment.
  646. """
  647. if len(waypoints) < 2:
  648. raise ValueError("build_waypoint_centers needs at least 2 waypoints.")
  649. centers: List[np.ndarray] = [waypoints[0]]
  650. for c0, c1 in zip(waypoints[:-1], waypoints[1:]):
  651. if centerline == 'curved':
  652. seg, _, _ = construct_centers_curved(atom_xyz, atom_r, c0, c1, step, curve_radius, curve_iters)
  653. else:
  654. seg = sample_straight_centers(c0, c1, step)
  655. centers.extend(seg[1:]) # skip seg[0]: it duplicates the previous segment's endpoint
  656. return centers
  659. def profile_along_centers(atom_xyz, atom_r, centers, eps, atom_meta, hydro_scale, electro_scale, occupancy_metric,
  660. true_area=True, probe=0.0, area_grid_step=0.25,
  661. area_min_radius=4.0, area_max_radius=20.0, area_radius_factor=3.0):
  662. rows = []
  663. s_acc = 0.0
  664. index = SpatialIndex(atom_xyz, atom_r)
  665. for i, c in enumerate(centers):
  666. if i > 0:
  667. s_acc += float(np.linalg.norm(centers[i] - centers[i - 1]))
  668. rmin, tags, hyd, elec = _evaluate_slice(
  669. atom_xyz, atom_r, atom_meta, c, eps, hydro_scale, electro_scale, index=index)
  670. rows.append({
  671. 's_A': float(s_acc), 'x': float(c[0]), 'y': float(c[1]), 'z': float(c[2]),
  672. 'radius_A': float(rmin), 'hydro_index': float(hyd), 'electro_index': float(elec),
  673. 'contributors': tags,
  674. })
  675. for i in range(len(rows)):
  676. if i == 0:
  677. t = np.array([rows[1]['x'] - rows[0]['x'], rows[1]['y'] - rows[0]['y'], rows[1]['z'] - rows[0]['z']], float)
  678. elif i == len(rows) - 1:
  679. t = np.array([rows[i]['x'] - rows[i - 1]['x'], rows[i]['y'] - rows[i - 1]['y'], rows[i]['z'] - rows[i - 1]['z']], float)
  680. else:
  681. t = np.array([rows[i + 1]['x'] - rows[i - 1]['x'], rows[i + 1]['y'] - rows[i - 1]['y'], rows[i + 1]['z'] - rows[i - 1]['z']], float)
  682. n = np.linalg.norm(t)
  683. t = np.array([0.0, 0.0, 1.0]) if n < 1e-9 else (t / n)
  684. rows[i]['tx'] = float(t[0])
  685. rows[i]['ty'] = float(t[1])
  686. rows[i]['tz'] = float(t[2])
  687. if true_area:
  688. # Cross-section area is measured perpendicular to the *local*
  689. # centerline direction, not a fixed global axis -- important once
  690. # the path is genuinely curved.
  691. v_local, w_local = orthonormal_basis_from_axis(t)
  692. c = np.array([rows[i]['x'], rows[i]['y'], rows[i]['z']], float)
  693. area_geom, area_access, half_extent = slice_area(
  694. index, c, v_local, w_local, probe=probe,
  695. min_half_extent=area_min_radius, max_half_extent=area_max_radius,
  696. radius_factor=area_radius_factor, grid_step=area_grid_step)
  697. rows[i]['area_A2'] = float(area_geom)
  698. rows[i]['area_access_A2'] = float(area_access)
  699. rows[i]['area_half_extent_A'] = float(half_extent)
  700. if occupancy_metric == 'hydro':
  701. rows[i]['occ_value'] = float(rows[i]['hydro_index'])
  702. elif occupancy_metric == 'electro':
  703. rows[i]['occ_value'] = float(rows[i]['electro_index'])
  704. else:
  705. rows[i]['occ_value'] = float(rows[i]['radius_A'])
  706. return rows
  709. def write_csv(path, rows):
  710. if not rows:
  711. return
  712. cols = list(rows[0].keys())
  713. if 'occ_value' not in cols:
  714. cols.append('occ_value')
  715. with open(path, 'w', newline='') as f:
  716. w = csv.DictWriter(f, fieldnames=cols)
  717. w.writeheader()
  718. for r in rows:
  719. w.writerow(r)
  722. def _format_pdb_atom_line(serial, name, resName, chainID, resSeq, x, y, z, occupancy, bfactor, element,
  723. altLoc=' ', iCode=' '):
  724. return (
  725. f"ATOM {serial:5d} {name:^4}{altLoc}{resName:>3} {chainID}{resSeq:>4}{iCode} "
  726. f"{x:8.3f}{y:8.3f}{z:8.3f}{occupancy:6.3f}{bfactor:6.3f} {element:>2} \r\n"
  727. )
  730. def write_mesh_pdb(path, rows, axis_u, rings=24):
  731. use_local = ('tx' in rows[0])
  732. if not use_local:
  733. u = axis_u / (np.linalg.norm(axis_u) + 1e-12)
  734. base_v, base_w = orthonormal_basis_from_axis(u)
  735. coords = []
  736. bvals = []
  737. occs = []
  738. for r in rows:
  739. c = np.array([r['x'], r['y'], r['z']], float)
  740. rad = max(MIN_PORE_RADIUS_A, float(r['radius_A']))
  741. occ = float(r.get('occ_value', 0.0))
  742. if use_local:
  743. u_loc = np.array([r['tx'], r['ty'], r['tz']], float)
  744. v, w = orthonormal_basis_from_axis(u_loc)
  745. else:
  746. v, w = base_v, base_w
  747. for k in range(rings):
  748. ang = 2 * np.pi * (k / rings)
  749. p = c + rad * (np.cos(ang) * v + np.sin(ang) * w)
  750. coords.append(p)
  751. bvals.append(rad)
  752. occs.append(occ)
  753. lines = []
  754. serial_start = 1
  755. for i, pnt in enumerate(coords, start=serial_start):
  756. x, y, z = pnt
  757. b = bvals[i - serial_start]
  758. occ = occs[i - serial_start]
  759. lines.append(_format_pdb_atom_line(i, 'C', 'ALA', 'M', 1, x, y, z, occ, b, 'C'))
  760. # ring & longitudinal CONECT like v8
  761. n = len(rows)
  762. R = rings
  764. def idx(step, k):
  765. return serial_start + step * R + k
  767. for step in range(n):
  768. for k in range(R):
  769. lines.append(f"CONECT{idx(step, k):5d}{idx(step, (k + 1) % R):5d}\r\n")
  770. for step in range(n - 1):
  771. for k in range(R):
  772. lines.append(f"CONECT{idx(step, k):5d}{idx(step + 1, k):5d}\r\n")
  773. with open(path, 'w', newline='') as f:
  774. f.writelines(lines)
  777. def write_centerline_pdb(path, rows):
  778. lines = []
  779. serial = 1
  780. for i, r in enumerate(rows, start=1):
  781. x, y, z = r['x'], r['y'], r['z']
  782. b = r['radius_A']
  783. occ = float(r.get('occ_value', 0.0))
  784. lines.append(f"HETATM{serial:5d} O PORE Z{i:4d} {x:8.3f}{y:8.3f}{z:8.3f}{occ:6.3f}{b:6.3f} O \r\n")
  785. serial += 1
  786. with open(path, 'w', newline='') as f:
  787. f.writelines(lines)
  790. # --- Extras (volumes, conductance, passability), matched to CLI v8 ---
  791. def _trap_volume_A3(rows, probe_A):
  792. if len(rows) < 2:
  793. return 0.0
  794. vol = 0.0
  795. for i in range(len(rows) - 1):
  796. ds = float(rows[i + 1]['s_A'] - rows[i]['s_A'])
  797. r0 = max(0.0, float(rows[i]['radius_A']) - probe_A)
  798. r1 = max(0.0, float(rows[i + 1]['radius_A']) - probe_A)
  799. vol += 0.5 * (math.pi * r0 * r0 + math.pi * r1 * r1) * ds
  800. return vol
  803. def _trap_true_volume_A3(rows, area_key: str) -> float:
  804. """Integrate a rasterized area column over s_A -- the true-cross-section
  805. counterpart to _trap_volume_A3's circular-radius approximation."""
  806. if len(rows) < 2 or area_key not in rows[0]:
  807. return 0.0
  808. s_vals = [r['s_A'] for r in rows]
  809. a_vals = [r[area_key] for r in rows]
  810. return _trapz(a_vals, s_vals)
  813. def _geometric_openness_index(rows, conductivity_S_per_m):
  814. """Unitless index of how geometrically open a pore is along its length --
  815. NOT a physical ionic conductance, and no longer reported in Siemens.
  817. Modeled loosely as the reciprocal of a series resistor-network integral
  818. (per-slice inverse cross-sectional area, plus an access-resistance-like
  819. end correction), using the parallel/series-resistor analogy for how a
  820. single bottleneck dominates a narrow pore's overall openness. This
  821. captures geometry only: it cannot account for ion occupancy,
  822. selectivity-filter state, channel gating, or electrostatics, so a
  823. wide-but-inactivated channel and a wide-and-conducting channel of
  824. identical shape score identically. ``conductivity_S_per_m`` is used only
  825. as an internal weighting term on the resistor-network terms above (it
  826. does not change the relative ranking of different pore geometries at a
  827. fixed conductivity value) -- it is not a claim about actual ionic
  828. conductance.
  829. """
  830. kappa = max(1e-12, float(conductivity_S_per_m))
  831. rho = 1.0 / kappa
  832. R_pore = 0.0
  833. blocked = False
  834. for i in range(len(rows) - 1):
  835. ds_m = float(rows[i + 1]['s_A'] - rows[i]['s_A']) * 1e-10
  836. r0 = max(1e-6, float(rows[i]['radius_A'])) * 1e-10
  837. r1 = max(1e-6, float(rows[i + 1]['radius_A'])) * 1e-10
  838. if r0 <= 0.0 or r1 <= 0.0:
  839. blocked = True
  840. A0 = math.pi * r0 * r0
  841. A1 = math.pi * r1 * r1
  842. R_pore += 0.5 * ((1.0 / max(A0, 1e-30)) + (1.0 / max(A1, 1e-30))) * ds_m * rho
  843. rmin_A = max(1e-6, float(min([r['radius_A'] for r in rows])))
  844. rmin_m = rmin_A * 1e-10
  845. R_access = (1.0 / (2.0 * kappa * max(rmin_m, 1e-12)))
  846. if blocked:
  847. return 0.0, float('inf'), R_access
  848. R_total = R_pore + R_access
  849. return (1.0 / R_total), R_pore, R_access
  852. def _passability_report(rows, pass_radii):
  853. pass_report = {}
  854. for sp, radA in pass_radii.items():
  855. radA = float(radA)
  856. spans = []
  857. start = None
  858. local_min = 1e9
  859. local_idx = None
  860. for i, row in enumerate(rows):
  861. rA = float(row['radius_A'])
  862. if rA < radA:
  863. if start is None:
  864. start = row['s_A']
  865. local_min = rA
  866. local_idx = i
  867. elif rA < local_min:
  868. local_min = rA
  869. local_idx = i
  870. elif start is not None:
  871. end = rows[i]['s_A']
  872. contrib = rows[local_idx]['contributors'] if local_idx is not None else ''
  873. spans.append({
  874. 'start_s_A': float(start), 'end_s_A': float(end),
  875. 'min_radius_A': float(local_min), 'min_contributors': contrib,
  876. })
  877. start = None
  878. local_min = 1e9
  879. local_idx = None
  880. if start is not None:
  881. end = rows[-1]['s_A']
  882. contrib = rows[local_idx]['contributors'] if local_idx is not None else ''
  883. spans.append({
  884. 'start_s_A': float(start), 'end_s_A': float(end),
  885. 'min_radius_A': float(local_min), 'min_contributors': contrib,
  886. })
  887. pass_report[sp] = {'is_passable': len(spans) == 0, 'blocked_spans': spans}
  888. return pass_report
  891. # --- CLI entry point ---
  892. def main(args=None):
  893. p = argparse.ArgumentParser(description="HOLE-like pore profile with straight/curved centerline options.")
  894. p.add_argument("pdb", help="Legacy PDB, or mmCIF/other gemmi-readable format (auto-converted)")
  895. p.add_argument("--top", default="")
  896. p.add_argument("--bottom", default="")
  897. p.add_argument("--mid", default="", help="Optional intermediate waypoint (same selection syntax as --top/--bottom). "
  898. "Routes the profile through this point between bottom and top -- useful "
  899. "for channels with a bend you already know the location of, without "
  900. "needing a full path search.")
  901. p.add_argument("--interactive", action="store_true")
  902. p.add_argument("--step", type=float, default=1.0)
  903. p.add_argument("--eps", type=float, default=0.25)
  904. p.add_argument("--noH", action="store_true")
  905. p.add_argument("--noHet", action="store_true")
  906. p.add_argument("--vdwjson", type=str, default="")
  907. p.add_argument("--rings", type=int, default=24)
  908. p.add_argument("--outprefix", type=str, default="holepy_out")
  909. p.add_argument("--probe", type=float, default=0.0)
  910. p.add_argument("--conductivity", type=float, default=1.5,
  911. help="Weights the internal resistor-network terms of the geometric openness index "
  912. "only; this is NOT a physical bulk-solution conductivity lookup and the index "
  913. "is not a reported conductance.")
  914. p.add_argument("--occupancy", choices=["hydro", "electro", "radii"], default="hydro")
  915. p.add_argument("--hydroscale", choices=["raw", "01"], default="raw")
  916. p.add_argument("--electroscale", choices=["raw", "01"], default="raw")
  917. p.add_argument("--passable_json", type=str, default="")
  918. p.add_argument("--centerline", choices=["straight", "curved"], default="straight")
  919. p.add_argument("--adaptive", action="store_true")
  920. p.add_argument("--slope_thresh", type=float, default=0.5)
  921. p.add_argument("--max_refine", type=int, default=3)
  922. p.add_argument("--curve_radius", type=float, default=2.0)
  923. p.add_argument("--curve_iters", type=int, default=3)
  924. p.add_argument("--no-true-area", dest="true_area", action="store_false",
  925. help="Skip rasterized true cross-sectional area/volume (faster; falls back to circular-radius approximation only)")
  926. p.add_argument("--area-grid-step", type=float, default=0.25, help="Grid spacing (A) for area rasterization")
  927. p.add_argument("--area-min-radius", type=float, default=4.0, help="Minimum half-extent (A) of the per-slice area search box")
  928. p.add_argument("--area-max-radius", type=float, default=20.0, help="Maximum half-extent (A) of the per-slice area search box")
  929. p.add_argument("--area-radius-factor", type=float, default=3.0, help="Area search box half-extent = this * local clearance, clamped to [min,max]")
  930. p.add_argument("--surface", action="store_true", help="Also write a density map (.mrc) for true isosurface rendering of the pore lumen")
  931. p.add_argument("--surface-voxel", type=float, default=0.5, help="Voxel spacing (A) for the density map")
  932. p.add_argument("--surface-margin", type=float, default=3.0, help="Extra padding (A) around the swept pore region in the density map")
  933. p.add_argument("--surface-radius-factor", type=float, default=1.3,
  934. help="Swept-tube half-width around the centerline for the density map = this * local clearance, "
  935. "clamped to [--area-min-radius, --area-max-radius]. Kept smaller/more conservative than "
  936. "--area-radius-factor to avoid picking up unrelated packing gaps ('swiss cheese') in the "
  937. "isosurface; increase only if the true surface looks clipped too close to the centerline.")
  938. p.set_defaults(true_area=True)
  939. ns = p.parse_args(args=args)
  941. try:
  942. pdb_path = ensure_legacy_pdb(Path(ns.pdb))
  943. except StructureTooLargeForLegacyPDB as e:
  944. print(f"ERROR: {e}", file=sys.stderr)
  945. return 4
  947. atoms = load_pdb_atoms(pdb_path, include_h=not ns.noH, include_hetatm=not ns.noHet)
  948. if ns.interactive or (not ns.top and not ns.bottom):
  949. print("Enter TOP:", file=sys.stderr)
  950. ns.top = input().strip()
  951. print("Enter BOTTOM:", file=sys.stderr)
  952. ns.bottom = input().strip()
  953. top_sel = parse_residue_tokens(ns.top)
  954. bot_sel = parse_residue_tokens(ns.bottom)
  955. ca_top = ca_positions_for(atoms, top_sel)
  956. ca_bot = ca_positions_for(atoms, bot_sel)
  957. if ca_top.size == 0 or ca_bot.size == 0:
  958. print("ERROR: Could not find CA atoms.", file=sys.stderr)
  959. return 3
  960. c_top = ca_top.mean(axis=0)
  961. c_bot = ca_bot.mean(axis=0)
  963. c_mid = None
  964. if ns.mid.strip():
  965. mid_sel = parse_residue_tokens(ns.mid)
  966. ca_mid = ca_positions_for(atoms, mid_sel)
  967. if ca_mid.size == 0:
  968. print("ERROR: Could not find CA atoms for --mid selection.", file=sys.stderr)
  969. return 3
  970. c_mid = ca_mid.mean(axis=0)
  972. custom_vdw = {}
  973. if ns.vdwjson:
  974. with open(ns.vdwjson, 'r') as jf:
  975. custom_vdw = normalize_vdw_keys(json.load(jf))
  977. coords = np.array([[a.x, a.y, a.z] for a in atoms], float)
  978. radii = np.array([vdw_radius(a.element, custom_vdw) for a in atoms], float)
  979. metas = [(a.chain, a.resname, a.resi, a.icode) for a in atoms]
  981. if c_mid is not None:
  982. waypoints = [c_bot, c_mid, c_top]
  983. centers = build_waypoint_centers(coords, radii, waypoints, ns.step, ns.centerline, ns.curve_radius, ns.curve_iters)
  984. rows = profile_along_centers(
  985. coords, radii, centers, ns.eps, metas, ns.hydroscale, ns.electroscale, ns.occupancy,
  986. true_area=ns.true_area, probe=ns.probe, area_grid_step=ns.area_grid_step,
  987. area_min_radius=ns.area_min_radius, area_max_radius=ns.area_max_radius,
  988. area_radius_factor=ns.area_radius_factor)
  989. if ns.adaptive:
  990. print("NOTE: --adaptive has no effect together with --mid (adaptive resampling only applies to a "
  991. "single straight segment); ignoring.", file=sys.stderr)
  992. elif ns.centerline == 'straight':
  993. rows, u, L = profile_along_axis(
  994. coords, radii, c_bot, c_top, ns.step, ns.eps, metas, adaptive=ns.adaptive,
  995. slope_thresh=ns.slope_thresh, max_refine=ns.max_refine,
  996. hydro_scale=ns.hydroscale, electro_scale=ns.electroscale, occupancy_metric=ns.occupancy,
  997. true_area=ns.true_area, probe=ns.probe, area_grid_step=ns.area_grid_step,
  998. area_min_radius=ns.area_min_radius, area_max_radius=ns.area_max_radius,
  999. area_radius_factor=ns.area_radius_factor)
  1000. else:
  1001. centers, u, L = construct_centers_curved(coords, radii, c_bot, c_top, ns.step, ns.curve_radius, ns.curve_iters)
  1002. rows = profile_along_centers(
  1003. coords, radii, centers, ns.eps, metas, ns.hydroscale, ns.electroscale, ns.occupancy,
  1004. true_area=ns.true_area, probe=ns.probe, area_grid_step=ns.area_grid_step,
  1005. area_min_radius=ns.area_min_radius, area_max_radius=ns.area_max_radius,
  1006. area_radius_factor=ns.area_radius_factor)
  1008. # volumes & conductance
  1009. volume_geom = _trap_volume_A3(rows, 0.0)
  1010. volume_access = _trap_volume_A3(rows, float(ns.probe))
  1011. volume_true_geom = _trap_true_volume_A3(rows, 'area_A2') if ns.true_area else None
  1012. volume_true_access = _trap_true_volume_A3(rows, 'area_access_A2') if ns.true_area else None
  1013. openness_index, R_pore, R_access = _geometric_openness_index(rows, ns.conductivity)
  1014. # passability reporting
  1015. pass_radii = {}
  1016. if ns.passable_json:
  1017. try:
  1018. with open(ns.passable_json, 'r') as pf:
  1019. pass_radii = json.load(pf)
  1020. except Exception:
  1021. pass
  1022. if not pass_radii:
  1023. pass_radii = {'water': 1.4, 'na': 1.02, 'k': 1.38, 'ca': 1.00}
  1024. pass_report = _passability_report(rows, pass_radii)
  1026. density_map_path = None
  1027. if ns.surface:
  1028. surf_index = SpatialIndex(coords, radii)
  1029. try:
  1030. values_zyx, origin, voxel = build_density_grid(
  1031. surf_index, rows, probe=ns.probe, voxel_size=ns.surface_voxel, margin=ns.surface_margin,
  1032. min_half_extent=ns.area_min_radius, max_half_extent=ns.area_max_radius,
  1033. radius_factor=ns.surface_radius_factor)
  1034. except ValueError as e:
  1035. print(f"ERROR: {e}", file=sys.stderr)
  1036. return 5
  1037. outprefix_tmp = Path(ns.outprefix)
  1038. if outprefix_tmp.suffix.lower() == '.csv':
  1039. outprefix_tmp = outprefix_tmp.with_suffix('')
  1040. density_map_path = outprefix_tmp.with_name(outprefix_tmp.stem + '_density.mrc')
  1041. write_density_map(density_map_path, origin, voxel, values_zyx)
  1043. outprefix = Path(ns.outprefix)
  1044. if outprefix.suffix.lower() == '.csv':
  1045. outprefix = outprefix.with_suffix('')
  1046. csv_path = outprefix.with_suffix('.csv')
  1047. pdb_center = outprefix.with_name(outprefix.stem + '_centerline.pdb')
  1048. pdb_mesh = outprefix.with_name(outprefix.stem + '_mesh.pdb')
  1049. summary_path = outprefix.with_name(outprefix.stem + '_summary.json')
  1051. with open(summary_path, 'w') as jf:
  1052. json.dump({
  1053. 'pdb': str(Path(ns.pdb).resolve()), 'top': ns.top, 'bottom': ns.bottom, 'mid': (ns.mid or None), 'centerline': ns.centerline,
  1054. 'adaptive': bool(ns.adaptive), 'slope_thresh': float(ns.slope_thresh), 'max_refine': int(ns.max_refine),
  1055. 'curve_radius_A': float(ns.curve_radius), 'curve_iters': int(ns.curve_iters),
  1056. 'step_A': float(ns.step), 'eps_A': float(ns.eps), 'probe_A': float(ns.probe),
  1057. 'length_A': float(rows[-1]['s_A'] if rows else 0.0),
  1058. 'volume_geometric_A3': float(volume_geom), 'volume_accessible_A3': float(volume_access),
  1059. 'volume_true_geometric_A3': (float(volume_true_geom) if volume_true_geom is not None else None),
  1060. 'volume_true_accessible_A3': (float(volume_true_access) if volume_true_access is not None else None),
  1061. 'true_area_enabled': bool(ns.true_area),
  1062. 'area_grid_step_A': float(ns.area_grid_step),
  1063. 'conductivity_S_per_m': float(ns.conductivity),
  1064. 'R_pore_ohm': float(R_pore), 'R_access_ohm': float(R_access), 'R_total_ohm': float(R_pore + R_access),
  1065. 'geometric_openness_index': float(openness_index),
  1066. 'passability': pass_report,
  1067. 'min_radius_A': float(min([r['radius_A'] for r in rows]) if rows else 0.0),
  1068. 'num_samples': len(rows),
  1069. 'occupancy_metric': ns.occupancy, 'hydroscale': ns.hydroscale, 'electroscale': ns.electroscale,
  1070. 'no_hydrogens': bool(ns.noH), 'ignore_hetatm': bool(ns.noHet),
  1071. 'density_map': (str(density_map_path) if density_map_path else None),
  1072. 'density_map_voxel_A': (float(ns.surface_voxel) if density_map_path else None),
  1073. }, jf, indent=2)
  1075. write_csv(csv_path, rows)
  1076. write_centerline_pdb(pdb_center, rows)
  1077. write_mesh_pdb(pdb_mesh, axis_u=(c_top - c_bot), rows=rows, rings=ns.rings)
  1078. print(f"Wrote: {csv_path}")
  1079. print(f"Wrote: {pdb_center}")
  1080. print(f"Wrote: {pdb_mesh}")
  1081. print(f"Wrote: {summary_path}")
  1082. if density_map_path:
  1083. print(f"Wrote: {density_map_path}")
  1084. print(f" -> In ChimeraX: open {density_map_path} then volume #<N> level 0")
  1085. print(f" (level 0 is the accessible pore boundary at probe radius {ns.probe:.2f} A)")
  1086. return 0
  1088. if __name__ == '__main__':
  1089. sys.exit(main())

pyHole.py at commit 890947b, under GPL-3.0 · at the source

Overview

Authors: Mariah R Baker1, Xiaoxuan Lin2, Guizhen Fan1, Ariel Martinez-Chavez1, Larry E Wagner2, Sundeep Malik2, Tyler Allison3, Briar Bell1, Alexander B Seryshev1, Julio Cordero-Morales1,3, Matthew L Baker1,3, David I Yule2, Irina I Serysheva1,3
ORCID iDs: Mariah R Baker, Xiaoxuan Lin, Guizhen Fan, Tyler Allison, Briar Bell, Julio Cordero-Morales, David I Yule, Irina I Serysheva
  1. Department of Biochemistry and Molecular Biology, McGovern Medical School, University of Texas Health Science Center, Houston, TX USA
  2. Department of Pharmacology and Physiology, University of Rochester, Rochester, NY USA
  3. MD Anderson Cancer Center, UTHealth Graduate School of Biomedical Sciences, University of Texas Health Science Center at Houston, Houston, TX USA
Institutions: The University of Texas MD Anderson Cancer Center (United States); The University of Texas Health Science Center at Houston (United States); University of Rochester (United States); University of Rochester Medicine (United States)
Journal: Nature communications, volume 17, issue 1, article 8929
Dates: received 18 December 2025; accepted 9 July 2026; published online 22 July 2026
Type: Research article · Language: English
License: CC BY-NC-ND
Identifiers: DOI 10.1038/s41467-026-75806-y · PMID 42637755 · PMCID PMC13503717 · OpenAlex W7170038317
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: histology / microscopy (modality), human (organism), cellular / molecular (subfield)
Methods: Spectral & time-frequency, fMRI & imaging
Keywords: Permeation and transport, Cryoelectron microscopy, Ion channels
MeSH: Cryoelectron Microscopy*, Inositol 1,4,5-Trisphosphate Receptors*, Adenosine Triphosphate, Animals, Binding Sites, Calcium, Humans, Inositol 1,4,5-Trisphosphate, Models, Molecular, Protein Conformation, Protein Isoforms (* major topic)
Topic: Calcium signaling and nucleotide metabolism (Physiology, Biochemistry, Genetics and Molecular Biology), according to OpenAlex
Funding: U.S. Department of Health &amp; Human Services | NIH | National Institute of General Medical Sciences (R21AR082833, R35GM153178, R35GM153178-01S1); U.S. Department of Health & Human Services | NIH | National Institute of General Medical Sciences (NIGMS) (R21AR082833, R35GM153178-01S1, R35GM153178); U.S. Department of Health & Human Services | NIH | NIH Office of the Director (OD) (S10OD032204); U.S. Department of Health &amp; Human Services | NIH | NIH Office of the Director (S10OD032204); NIH HHS (S10 OD032204); NIAMS NIH HHS (R21 AR082833); NIGMS NIH HHS (R35 GM153178)
Citations: not cited yet (Europe PMC); 77 references in the paper

Abstract

The abstract is not reproduced here: the paper's license (CC BY-NC-ND) does not allow it. Read it in the paper, at the publisher or on Europe PMC.

Repository

Its files are read in the Code ↔ Paper reader above, with 2 matches between paragraphs and lines of code.

mlbaker-uth/pyHole

License: GPL-3.0
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 890947b36475fa200411536bb929bdd8dac7cbea, 23 September 2026
Languages: Python (15)
Size: 31 files, 15 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: README, license file, environment (requirements.txt)
Not found: CITATION.cff, tests, continuous integration, documentation
Tools: NumPy (6 files), SciPy (3 files), Matplotlib (1 file), pandas (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
17 files
At the source: github.com/mlbaker-uth/pyHole

Code availability statement

The paper has a code availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

Read it in the paper: doi.org/10.1038/s41467-026-75806-y.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 15 scripts, each with its path and the digest of its content;
  • 2 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

Code and data availability statement

The paper has a code and data availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

Read it in the paper: doi.org/10.1038/s41467-026-75806-y.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 13 authors, 3 keywords, 11 MeSH terms, 7 funders, 76 references.

Cite

This paper

Baker, M. R., Lin, X., Fan, G., Martinez-Chavez, A., Wagner, L. E., Malik, S., Allison, T., Bell, B., Seryshev, A. B., Cordero-Morales, J., Baker, M. L., Yule, D. I., & Serysheva, I. I. (2026). Cryo-EM insights into isoform-specific properties of the IP&lt;sub&gt;3&lt;/sub&gt;R2 channel. Nature communications, 17(1), 8929. https://doi.org/10.1038/s41467-026-75806-y

BibTeX

@article{baker2026cryo,
author = {Baker, Mariah R and Lin, Xiaoxuan and Fan, Guizhen and Martinez-Chavez, Ariel and Wagner, Larry E and Malik, Sundeep and Allison, Tyler and Bell, Briar and Seryshev, Alexander B and Cordero-Morales, Julio and Baker, Matthew L and Yule, David I and Serysheva, Irina I},
title = {{Cryo-EM insights into isoform-specific properties of the IP\&lt;sub\&gt;3\&lt;/sub\&gt;R2 channel}},
journal = {Nature communications},
year = {2026},
month = jul,
volume = {17},
number = {1},
pages = {8929},
publisher = {Nature Publishing Group},
issn = {2041-1723},
doi = {10.1038/s41467-026-75806-y},
url = {https://doi.org/10.1038/s41467-026-75806-y},
pmid = {42637755},
pmcid = {PMC13503717}
}

RIS

TY - JOUR
AU - Baker, Mariah R
AU - Lin, Xiaoxuan
AU - Fan, Guizhen
AU - Martinez-Chavez, Ariel
AU - Wagner, Larry E
AU - Malik, Sundeep
AU - Allison, Tyler
AU - Bell, Briar
AU - Seryshev, Alexander B
AU - Cordero-Morales, Julio
AU - Baker, Matthew L
AU - Yule, David I
AU - Serysheva, Irina I
TI - Cryo-EM insights into isoform-specific properties of the IP&lt;sub&gt;3&lt;/sub&gt;R2 channel
T2 - Nature communications
J2 - Nat Commun
PY - 2026
DA - 2026/07/22
VL - 17
IS - 1
SP - 8929
SN - 2041-1723
PB - Nature Publishing Group
DO - 10.1038/s41467-026-75806-y
UR - https://doi.org/10.1038/s41467-026-75806-y
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41467-026-75806-y",
"type": "article-journal",
"title": "Cryo-EM insights into isoform-specific properties of the IP&lt;sub&gt;3&lt;/sub&gt;R2 channel",
"container-title": "Nature communications",
"author": [
{
"family": "Baker",
"given": "Mariah R"
},
{
"family": "Lin",
"given": "Xiaoxuan"
},
{
"family": "Fan",
"given": "Guizhen"
},
{
"family": "Martinez-Chavez",
"given": "Ariel"
},
{
"family": "Wagner",
"given": "Larry E"
},
{
"family": "Malik",
"given": "Sundeep"
},
{
"family": "Allison",
"given": "Tyler"
},
{
"family": "Bell",
"given": "Briar"
},
{
"family": "Seryshev",
"given": "Alexander B"
},
{
"family": "Cordero-Morales",
"given": "Julio"
},
{
"family": "Baker",
"given": "Matthew L"
},
{
"family": "Yule",
"given": "David I"
},
{
"family": "Serysheva",
"given": "Irina I"
}
],
"container-title-short": "Nat Commun",
"volume": "17",
"issue": "1",
"page": "8929",
"DOI": "10.1038/s41467-026-75806-y",
"PMID": "42637755",
"PMCID": "PMC13503717",
"ISSN": "2041-1723",
"publisher": "Nature Publishing Group",
"URL": "https://doi.org/10.1038/s41467-026-75806-y",
"language": "en",
"issued": {
"date-parts": [
[
2026,
7,
22
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s41467-026-75444-4 [code]
Structural insights enable drug discovery for the neuronal NBCn2 carbonate transporter.
Journal: Nature communications
In common: pandas, SciPy, Matplotlib, 1 other tool, ebi.ac.uk/pdbe/entry, histology / microscopy, cellular / molecular, 5 references
[2] doi:10.1038/s41467-026-74087-9 [code]
Cryo-EM structures of heteromeric Kir4.1/5.1 channel suggest mechanisms of inward rectification and channel blockage.
Journal: Nature communications
In common: ebi.ac.uk/pdbe/entry, histology / microscopy, cellular / molecular, 6 references
[3] doi:10.1038/s41594-026-01866-9 [code]
Structural and mechanistic insights into gating and allosteric modulation of GluN1-GluN3A NMDA receptors.
Journal: Nature structural & molecular biology
In common: SciPy, Matplotlib, NumPy, ebi.ac.uk/pdbe/entry, histology / microscopy, cellular / molecular, 4 references
[4] doi:10.1038/s41594-026-01845-0
Conformational plasticity of human acid-sensing ion channel 1a.
Journal: Nature structural & molecular biology
In common: ebi.ac.uk/pdbe/entry, histology / microscopy, cellular / molecular, 4 references
[5] doi:10.1371/journal.pbio.3003777
Structure of the human P2X3 receptor reveals the basis for subtype-selective inhibition by sivopixant.
Journal: PLoS biology
In common: ebi.ac.uk/pdbe/entry, histology / microscopy, cellular / molecular, 4 references
[6] doi:10.1038/s41467-026-74814-2
Structural mechanism of Necrocide 1 activation of human TRPM4 that triggers necrosis by sodium overload.
Journal: Nature communications
In common: ebi.ac.uk/pdbe/entry, histology / microscopy, cellular / molecular, 4 references
[7] doi:10.1038/s41467-026-70575-0
Structure of a pH-sensitive pentameric ligand-gated ion channel from the Sarcoptes scabies mite.
Journal: Nature communications
In common: ebi.ac.uk/pdbe/entry, histology / microscopy, cellular / molecular, 4 references
[8] doi:10.1038/s41594-026-01787-7 [code]
Microtubules in the axon are GDP bound but adopt a stable GTP-like expanded state.
Journal: Nature structural & molecular biology
In common: Matplotlib, NumPy, ebi.ac.uk/pdbe/entry, histology / microscopy, cellular / molecular, 2 references
[9] doi:10.1038/s41467-026-75564-x
Cooperative mechanism of neurotransmitter recognition and transport by the human vesicular polyamine transporter.
Journal: Nature communications
In common: ebi.ac.uk/pdbe/entry, histology / microscopy, cellular / molecular, 3 references
[10] doi:10.1038/s41467-026-75877-x
Structure of NHE6 and its lipid-mediated interactions regulating endosomal pH.
Journal: Nature communications
In common: ebi.ac.uk/pdbe/entry, histology / microscopy, cellular / molecular, 3 references

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.