molframe_spatial/
base_pair_geometry.rs

1//! Dictionary-free canonical pairing geometry, O(residues + local candidates).
2
3use crate::{CellList, SpatialError};
4use molframe_chem::base_pair::{self as edges, CanonicalBase, WatsonCrickGeometry};
5use molframe_core::{ResidueIndex, Structure};
6use num_traits::ToPrimitive;
7
8/// Canonical edge contacts and coplanar sugar geometry supporting a pair.
9/// This is geometric evidence; it does not report perceived hydrogen bonds.
10#[derive(Clone, Copy, Debug, PartialEq)]
11pub struct GeometricBasePair {
12    /// Lower-indexed nucleotide residue.
13    pub first: ResidueIndex,
14    /// Higher-indexed nucleotide residue.
15    pub second: ResidueIndex,
16    /// Number of canonical edge contacts within the requested distance.
17    pub contact_count: usize,
18    /// Shortest canonical edge contact, in ångström.
19    pub closest_distance: f32,
20}
21
22/// Finds standard A–T/U and G–C pairs without a component dictionary.
23///
24/// Requires positioned C1', C2, C4 and C6 atoms, reference-plane agreement,
25/// and at least two complementary edge contacts. Unknown or modified component
26/// names are skipped; dictionary-backed chemistry is available through
27/// the `molframe::analysis` domain. Positions are read from the first model.
28///
29/// # Errors
30/// Returns a spatial error for non-finite or non-positive distance limits.
31pub fn standard_base_pair_geometry(
32    structure: &Structure,
33    geometry: WatsonCrickGeometry,
34    maximum_contact_distance: f32,
35) -> Result<Vec<GeometricBasePair>, SpatialError> {
36    if !maximum_contact_distance.is_finite()
37        || maximum_contact_distance <= 0.0
38        || !geometry.minimum_c1_distance.is_finite()
39        || geometry.minimum_c1_distance < 0.0
40        || !geometry.maximum_c1_distance.is_finite()
41        || geometry.maximum_c1_distance <= geometry.minimum_c1_distance
42        || !geometry.maximum_plane_angle_degrees.is_finite()
43        || !(0.0..=90.0).contains(&geometry.maximum_plane_angle_degrees)
44    {
45        return Err(SpatialError::InvalidCutoff);
46    }
47    let mut residues = Vec::new();
48    let mut positions = Vec::new();
49    for residue in structure.data().residues() {
50        let Some(base) = residue.name().and_then(standard_base) else {
51            continue;
52        };
53        let Some(position) = residue
54            .atom("C1'")
55            .and_then(molframe_core::structure::AtomRef::position)
56        else {
57            continue;
58        };
59        residues.push((residue, base));
60        positions.push(position);
61    }
62    let rows: Vec<_> = (0..positions.len())
63        .map(u32::try_from)
64        .collect::<Result<_, _>>()
65        .map_err(|_| SpatialError::InvalidCutoff)?;
66    let Some(cutoff) = geometry.maximum_c1_distance.to_f32() else {
67        return Err(SpatialError::InvalidCutoff);
68    };
69    let index = CellList::build(&positions, &rows, cutoff, None)?;
70    let mut pairs = Vec::new();
71    for (row, &(left, left_base)) in residues.iter().enumerate() {
72        index.for_each_neighbor(positions[row], cutoff, |other, _| {
73            let other = other as usize;
74            if other <= row {
75                return;
76            }
77            let (right, right_base) = residues[other];
78            if !left_base.complementary(right_base)
79                || !edges::has_watson_crick_geometry(left, right, geometry)
80            {
81                return;
82            }
83            let mut count = 0;
84            let mut closest = f32::INFINITY;
85            for first in ["N1", "N2", "N3", "N4", "N6", "O2", "O4", "O6"] {
86                for second in ["N1", "N2", "N3", "N4", "N6", "O2", "O4", "O6"] {
87                    if !edges::watson_crick_contact(left_base, first, right_base, second) {
88                        continue;
89                    }
90                    let (Some(a), Some(b)) = (
91                        left.atom(first)
92                            .and_then(molframe_core::structure::AtomRef::position),
93                        right
94                            .atom(second)
95                            .and_then(molframe_core::structure::AtomRef::position),
96                    ) else {
97                        continue;
98                    };
99                    let squared: f32 = a.into_iter().zip(b).map(|(a, b)| (a - b).powi(2)).sum();
100                    if squared <= maximum_contact_distance.powi(2) {
101                        count += 1;
102                        closest = closest.min(squared.sqrt());
103                    }
104                }
105            }
106            if count >= 2 {
107                pairs.push(GeometricBasePair {
108                    first: left.index(),
109                    second: right.index(),
110                    contact_count: count,
111                    closest_distance: closest,
112                });
113            }
114        })?;
115    }
116    pairs.sort_by_key(|pair| (pair.first, pair.second));
117    Ok(pairs)
118}
119
120fn standard_base(name: &str) -> Option<CanonicalBase> {
121    match name {
122        "A" | "DA" => Some(CanonicalBase::Adenine),
123        "C" | "DC" => Some(CanonicalBase::Cytosine),
124        "G" | "DG" => Some(CanonicalBase::Guanine),
125        "T" | "DT" => Some(CanonicalBase::Thymine),
126        "U" | "DU" => Some(CanonicalBase::Uracil),
127        _ => None,
128    }
129}
130
131#[cfg(test)]
132#[path = "base_pair_geometry_tests.rs"]
133mod tests;
molframe_spatial/base_pair_geometry.rs

molframe_spatial/
base_pair_geometry.rs
