topo.analysis package

Submodules

topo.analysis.native_contacts module

Per-domain and per-interface fraction of native contacts (Q).

Computes the fraction of native contacts Q for every domain (intra-domain) and every domain-domain interface (inter-domain) declared in a domain.yaml, for each frame of a coarse-grained (CG, Cα-only) trajectory.

Method

  1. Native contacts are defined ONCE from the all-atom reference structure: residues i, j form a native contact when they are more than local_separation residues apart in sequence (abs(i - j) > 3 by default) AND at least one pair of heavy atoms (one per residue) is within cutoff (4.5 Å by default). The reference Cα–Cα distance of each contact is stored.

  2. Contacts are grouped by domain using the residue sets in domain.yaml: - intra-domain : both residues in the same domain - interface : one residue in each of two domains

  3. For every CG frame, a native contact (i, j) is counted as formed when the instantaneous Cα–Cα distance satisfies d_CG(i, j) <= tolerance * d_native(i, j) (tolerance = 1.2). Q is the fraction of that group’s native contacts formed.

Use as a library (e.g. from the nscale optimizer) by importing the functions, or as a CLI:

python -m topo.analysis.native_contacts -d domain.yaml -r ref.pdb

-p cg.psf -f cg.dcd -o Q.csv

Residue numbering: domain.yaml uses 1-based residue numbers that run 1..n_residues in CG bead order; internally these map to 0-based indices (n - 1), which is also the Cα/bead order of the reference and CG models.

topo.analysis.native_contacts.build_native_contacts(group1, group2, heavy_positions, heavy_res, ca_positions, cutoff, local_separation, disorder=None)[source]

Build the native-contact list for one domain or interface.

Parameters:

  • group1 (set[int] or None) – 0-based residue-index sets: - group1 is None, group2 is None -> all residues (whole protein) - group1 set, group2 is None -> intra-domain (both in group1) - group1 set, group2 set -> interface (one in each)

  • group2 (set[int] or None) – 0-based residue-index sets: - group1 is None, group2 is None -> all residues (whole protein) - group1 set, group2 is None -> intra-domain (both in group1) - group1 set, group2 set -> interface (one in each)

  • heavy_positions ((A, 3) array) – Coordinates of all reference heavy (non-H) atoms.

  • heavy_res ((A,) int array) – 0-based residue index of each heavy atom.

  • ca_positions ((N, 3) array) – Reference Cα coordinates in 0-based residue order.

  • cutoff (float) – Heavy-atom distance defining a contact (Å).

  • local_separation (int) – Minimum sequence separation (abs(i - j) > local_separation).

  • disorder (set[int] or None, optional) – 0-based indices of disordered (IDR) residues. Any native contact with either residue in this set is dropped, matching the energy path where IDR pairs carry no native contact (spec §2.7). None (default) keeps the pre-IDR behavior, so other callers are byte-identical.

Returns:

  • pairs ((M, 2) int array) – Native-contact residue index pairs (i < j), 0-based.

  • native_dist ((M,) float array) – Reference Cα–Cα distance for each contact.

topo.analysis.native_contacts.fraction_native_contacts(ca_positions, pairs, native_dist, tolerance)[source]

Q for one CG frame: fraction of native contacts with d_CG <= tolerance * d_native. Returns NaN if the group has no contacts.

topo.analysis.native_contacts.load_domains(domain_yaml, n_res)[source]

Read domains from domain.yaml as 0-based residue-index sets.

An optional disordered: section marks IDR residues (spec §2.7): those residues are subtracted from every domain (effective membership = domain - disordered), mirroring the energy path where “disorder wins”, and returned as a set so the caller can also drop any IDR-involving native contact from Q. intra_domains is optional (spec §2.8): a fully-IDP domain_def is just disordered: + n_residues and yields an empty domains dict.

Returns:

  • domains (dict[str, np.ndarray]) – Domain name -> sorted array of 0-based residue indices (IDR residues removed).

  • n_residues (int) – n_residues declared in the file (validated against the structure).

  • disorder (set[int]) – 0-based indices of disordered residues (empty when no disordered: section).

topo.analysis.native_contacts.load_universe(*args, **kwargs)[source]

mda.Universe with the harmless placeholder-CRYST1 warning suppressed.

Some PDBs carry a dummy CRYST1  1 1 1 record; MDAnalysis then warns “1 A^3 CRYST1 record, this is usually a placeholder. Unit cell dimensions will be set to None.” That is expected for these non-periodic CG inputs, so we silence only that one UserWarning and let any others through.

topo.analysis.native_contacts.main()[source]

Score per-domain and per-interface native contacts for a CG trajectory.

Builds the native-contact lists once from the all-atom reference (whole protein, each domain, and each domain pair interface), then computes the fraction of native contacts (Q) for every frame of the CG DCD trajectory and writes the per-frame results to the output CSV.

Raises:

  • FileNotFoundError – If any of the reference, PSF, DCD or domain input files is missing.

  • ValueError – If the domain definition or CG bead count does not match the reference protein’s residue count.

topo.analysis.native_contacts.parse_args()[source]

Parse command-line arguments for the native-contacts (Q) driver.

Returns:

Parsed arguments: domain, reference, psf, dcd, output paths plus the cutoff, local_separation, tolerance, start_frame and end_frame scoring options.

Return type:

argparse.Namespace

topo.analysis.native_contacts.reference_residue_geometry(ref_universe)[source]

Return reference Cα positions (0-based residue order) and a mapping from each heavy atom to its 0-based protein-residue index.

Building an explicit universe-resindex -> 0-based position map makes the indexing robust to hetero residues or non-sequential numbering.

Module contents

Post-simulation trajectory analysis for TOPO CG models.

Currently provides native-contact order parameters (Q), per domain and per interface; future order parameters (dRMS, RMSD, inter-domain distances) belong here too.

topo.analysis.build_native_contacts(group1, group2, heavy_positions, heavy_res, ca_positions, cutoff, local_separation, disorder=None)[source]

Build the native-contact list for one domain or interface.

Parameters:

  • group1 (set[int] or None) – 0-based residue-index sets: - group1 is None, group2 is None -> all residues (whole protein) - group1 set, group2 is None -> intra-domain (both in group1) - group1 set, group2 set -> interface (one in each)

  • group2 (set[int] or None) – 0-based residue-index sets: - group1 is None, group2 is None -> all residues (whole protein) - group1 set, group2 is None -> intra-domain (both in group1) - group1 set, group2 set -> interface (one in each)

  • heavy_positions ((A, 3) array) – Coordinates of all reference heavy (non-H) atoms.

  • heavy_res ((A,) int array) – 0-based residue index of each heavy atom.

  • ca_positions ((N, 3) array) – Reference Cα coordinates in 0-based residue order.

  • cutoff (float) – Heavy-atom distance defining a contact (Å).

  • local_separation (int) – Minimum sequence separation (abs(i - j) > local_separation).

  • disorder (set[int] or None, optional) – 0-based indices of disordered (IDR) residues. Any native contact with either residue in this set is dropped, matching the energy path where IDR pairs carry no native contact (spec §2.7). None (default) keeps the pre-IDR behavior, so other callers are byte-identical.

Returns:

  • pairs ((M, 2) int array) – Native-contact residue index pairs (i < j), 0-based.

  • native_dist ((M,) float array) – Reference Cα–Cα distance for each contact.

topo.analysis.chirality_agreement(chi_frames, chi_ref)[source]

Fraction agreement K per frame: fraction of windows whose chirality sign matches the native reference. K -> 1 native-handed, K -> 0 mirror.

chi_frames may be (n_windows,) for a single frame or (n_frames, n_windows). Returns a (n_frames,) array.

topo.analysis.classify_mirror(Q, K, rmsd_native, rmsd_reflected, q_fold=0.5, k_thresh=0.3, rmsd_ratio_thresh=0.9, eps=1e-12)[source]

Per-frame mirror classification.

Returns:

  • rmsd_ratio (ndarray) – rmsd_reflected / rmsd_native (denominator floored by eps).

  • is_mirror (bool ndarray) – (Q > q_fold) AND (K < k_thresh) AND (rmsd_ratio < rmsd_ratio_thresh).

topo.analysis.fraction_native_contacts(ca_positions, pairs, native_dist, tolerance)[source]

Q for one CG frame: fraction of native contacts with d_CG <= tolerance * d_native. Returns NaN if the group has no contacts.

topo.analysis.kabsch_proper(mobile, target)[source]

Rotate centered mobile onto centered target (proper rotation only).

Both inputs must already be centered. The diag(1, 1, d) determinant correction (d = sign(det(V @ Wt))) forces det(rotation) = +1 so the fit can never silently absorb a reflection – the crux of mirror detection.

Returns the rotated mobile coordinates (not translated).

topo.analysis.load_domains(domain_yaml, n_res)[source]

Read domains from domain.yaml as 0-based residue-index sets.

An optional disordered: section marks IDR residues (spec §2.7): those residues are subtracted from every domain (effective membership = domain - disordered), mirroring the energy path where “disorder wins”, and returned as a set so the caller can also drop any IDR-involving native contact from Q. intra_domains is optional (spec §2.8): a fully-IDP domain_def is just disordered: + n_residues and yields an empty domains dict.

Returns:

  • domains (dict[str, np.ndarray]) – Domain name -> sorted array of 0-based residue indices (IDR residues removed).

  • n_residues (int) – n_residues declared in the file (validated against the structure).

  • disorder (set[int]) – 0-based indices of disordered residues (empty when no disordered: section).

topo.analysis.local_chirality(points, eps=1e-08)[source]

Normalized local chirality along an ordered point sequence.

For bond vectors v_k = P[k+1] - P[k] and each window i:

chi_i = ((v_i x v_{i+1}) . v_{i+2}) / (|v_i| |v_{i+1}| |v_{i+2}|)

a signed pseudo-scalar in [-1, 1] (handedness). It flips sign under a reflection and is invariant to proper rotation and translation.

Parameters:

  • points ((M, 3) or (n_frames, M, 3) array) – Ordered points (the SS segment-endpoint polyline). A single (M, 3) frame is promoted internally and returned as (M - 3,).

  • eps (float) – Added to the denominator to guard zero-length bond vectors (coincident endpoints).

Returns:

chi

Return type:

(M - 3,) or (n_frames, M - 3) array

topo.analysis.mirror_rmsd(ca_frame, ca_native, reflect_axis='x')[source]

Per-frame RMSD to the native and to the reflected-native reference.

Parameters:

  • ca_frame ((N, 3) array) – Cα coordinates of the frame.

  • ca_native ((N, 3) array) – Cα coordinates of the native reference (any centering; centered here).

  • reflect_axis ({"x", "y", "z"}) – Axis negated to build the reflected reference. Any single axis works – the proper-rotation fit re-orients freely, so the choice does not affect the result.

Returns:

(rmsd_native, rmsd_reflected) – rmsd_reflected < rmsd_native indicates a mirror-like fold.

Return type:

tuple of float

topo.analysis.reference_residue_geometry(ref_universe)[source]

Return reference Cα positions (0-based residue order) and a mapping from each heavy atom to its 0-based protein-residue index.

Building an explicit universe-resindex -> 0-based position map makes the indexing robust to hetero residues or non-sequential numbering.

topo.analysis.segment_endpoints_from_file(segments_file)[source]

0-based Cα bead indices from an explicit <label> <start> <end> segment file (1-based residue numbers, one element per line), flattened and sorted.

Residue numbers are taken to run 1..n_residues in CG bead order (single-chain convention, as in the standalone prototypes); use the STRIDE path for multi-chain structures.

topo.analysis.segment_endpoints_from_stride(stride_output_file, key_to_index, min_len=4)[source]

0-based Cα bead indices of the head & tail of each STRIDE helix/strand element (length >= min_len), flattened across elements and sorted -> the ordered point sequence for the chirality polyline.

Reuses topo’s (chain, resid) -> 0-based index map so the endpoints line up with the CG bead order used everywhere else (and multi-chain structures resolve).

STRIDE LOC record (fixed columns), e.g.:

LOC  AlphaHelix   MET     1 A      MET     13 A
LOC  Strand       SER    20 A      PRO     27 A