Merge pull request #16 from TKanX/bugfix/15-nucleic-acid-5-terminal-repair-and-hydrogenation-mishandles-phosphate-groups

fix(ops): Add Support for 5'-Terminal Nucleic Acid Phosphates

Author: TKanX

ARCHITECTURE.md

@@ -78,7 +78,7 @@ flowchart TD
 - **FetchTemplate** – pulls `TemplateView` for each residue name, failing early if missing.
 - **AlignPairs** – pairs existing heavy atoms with template coordinates for alignment anchors.
 - **SVD** – computes rotation/translation via Kabsch/SVD, handling single/dual point fallbacks.
-- **SynthesizeAtoms** – recreates missing heavy atoms (and terminal OXT when needed) using transformed template positions.
+- **SynthesizeAtoms** – recreates missing heavy atoms (terminal OXT for peptides, OP3 for 5'-phosphorylated nucleic acids) using transformed template positions or tetrahedral geometry.
 - **Cleanup** – removes atoms not present in the template to ensure canonical composition.
 
 ### 3.3 Hydro Pipeline (`ops::hydro`)
@@ -98,7 +98,7 @@ flowchart TD
 - **Protonation** – applies pH-driven heuristics plus HIS strategy selection to decide residue names.
 - **StripOldH** – removes existing hydrogens if `remove_existing_h` is enabled.
 - **GeometryBuild** – reconstructs hydrogens using template anchors, calling `reconstruct_geometry` and `calculate_transform` internally.
-- **TerminalAdjust** – adds terminal hydrogens (N-term, C-term, nucleic 5'/3') respecting protonation states.
+- **TerminalAdjust** – adds terminal hydrogens (N-term H1/H2/H3, C-term HOXT, nucleic 5' HO5' or pH-dependent phosphate HOP3, nucleic 3' HO3') respecting protonation states.
 - **EndHydro** – structure now contains geometrically sound hydrogens.
 
 ### 3.4 Solvate Pipeline (`ops::solvate`)
@@ -136,7 +136,7 @@ flowchart TD
 - **PrepareOffsets** – computes atom-index offsets for every residue to map local indexes to global topology indexes.
 - **IntraStandard** – connects intra-residue bonds per template definitions, including hydrogen anchors.
 - **HeteroTemplates** – injects bonds for hetero residues via user-provided `Template`s or errors if absent.
-- **TerminalBonds** – adds special-case bonds for terminal atoms (H1/H2/H3, HOXT, HO5', HO3').
+- **TerminalBonds** – adds special-case bonds for terminal atoms (H1/H2/H3, HOXT for peptides; P–OP3, OP3–HOP3, O5'–HO5', O3'–HO3' for nucleic acids).
 - **InterResidue** – detects peptide and nucleic linkages by measuring atom distances against cutoffs.
 - **DisulfideScan** – adds bonds between cystine sulfurs within the disulfide cutoff.
 - **EmitTopology** – produces the final `Topology` pairing the structure with collected bonds.
@@ -179,7 +179,7 @@ flowchart TD
 - **Anchor selection** – each template hydrogen lists one or more anchor atoms; missing anchors trigger `IncompleteResidueForHydro` errors to avoid guesswork.
 - **Rigid transform** – `reconstruct_geometry` retrieves the residue-specific transform (rotation + translation) derived from current heavy atoms and applies it to the template hydrogen coordinate.
 - **Randomization** – none is applied for standard hydrogens, ensuring deterministic placement; terminals use evenly spaced tetrahedral vectors sorted by dot product to preserve orientation.
-- **Terminal logic** – N-termini place up to three hydrogens arranged around the N–CA axis, while C-termini and nucleic terminals enforce the specific HOXT/HO5'/HO3' connectivity.
+- **Terminal logic** – N-termini place up to three hydrogens arranged around the N–CA axis; C-termini add HOXT to OXT; nucleic 5'-terminals either add HO5' (no phosphate) or pH-dependent HOP3 (with phosphate, below pKₐ₂ ≈ 6.5); nucleic 3'-terminals always add HO3'.
 
 ### 4.3 Ion Replacement and Degradation Handling
 

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "bio-forge"
-version = "0.2.1"
+version = "0.2.2"
 authors = ["Tony Kan <tianchengkan@gmail.com>"]
 license = "MIT"
 repository = "https://github.com/TKanX/bio-forge"

MANUAL.md

@@ -72,7 +72,7 @@ Options:
 bioforge repair -i cleaned.pdb -o repaired.pdb
 ```
 
-- Aligns each standard residue to its template and fills in missing heavy atoms, including peptide termini.
+- Aligns each standard residue to its template and fills in missing heavy atoms, including peptide termini (OXT) and nucleic acid 5'-terminal phosphate (OP3).
 - Ideal immediately after `clean` to ensure the structure is chemically complete before protonation.
 
 ### `hydro` – Add hydrogens with titration awareness

README.md

@@ -14,7 +14,7 @@
 
 1. **Load** – `io::read_pdb_structure` or `io::read_mmcif_structure` parses coordinates with help from `IoContext` alias resolution.
 2. **Clean** – `ops::clean_structure` removes waters, ions, hetero residues, or arbitrary residue names via `CleanConfig`.
-3. **Repair** – `ops::repair_structure` realigns residues to their templates and rebuilds missing heavy atoms (including OXT on C-termini).
+3. **Repair** – `ops::repair_structure` realigns residues to their templates and rebuilds missing heavy atoms (including OXT on C-termini and OP3 on 5'-phosphorylated nucleic acids).
 4. **Hydrogenate** – `ops::add_hydrogens` infers protonation states (configurable pH and histidine strategy) and reconstructs hydrogens from template anchors.
 5. **Solvate/Ionize** – `ops::solvate_structure` creates a periodic box, packs water on a configurable lattice, and swaps molecules for ions to satisfy a target charge.
 6. **Topology** – `ops::TopologyBuilder` replays template bond definitions, peptide-link detection, nucleic backbone connectivity, and disulfide heuristics to emit a `Topology` object.
@@ -44,7 +44,7 @@ BioForge is also available as a library crate. Add it to your `Cargo.toml` depen
 
 ```toml
 [dependencies]
-bio-forge = "0.2.1"
+bio-forge = "0.2.2"
 ```
 
 #### Example: Preparing a PDB Structure

src/ops/hydro.rs

@@ -22,6 +22,8 @@ const DISULFIDE_SG_THRESHOLD: f64 = 2.2;
 const N_TERM_PKA: f64 = 8.0;
 /// Henderson–Hasselbalch breakpoint for protonated C-termini.
 const C_TERM_PKA: f64 = 3.1;
+/// Henderson–Hasselbalch breakpoint for the second dissociation of terminal phosphate.
+const PHOSPHATE_PKA2: f64 = 6.5;
 
 /// Parameters controlling hydrogen addition behavior.
 ///
@@ -370,7 +372,9 @@ fn construct_hydrogens_for_residue(
             construct_3_prime_hydrogen(residue)?;
         }
         ResiduePosition::FivePrime if residue.standard_name.is_some_and(|s| s.is_nucleic()) => {
-            if !residue.has_atom("P") && residue.has_atom("O5'") {
+            if residue.has_atom("P") {
+                construct_5_prime_phosphate_hydrogens(residue, config)?;
+            } else if residue.has_atom("O5'") {
                 construct_5_prime_hydrogen(residue)?;
             }
         }
@@ -654,6 +658,56 @@ fn construct_5_prime_hydrogen(residue: &mut Residue) -> Result<(), Error> {
     Ok(())
 }
 
+/// Adds hydrogens to 5'-terminal phosphate groups based on pH.
+///
+/// At physiological pH (≥6.5), the terminal phosphate carries two negative charges and
+/// requires no protons. Below this threshold, one proton is added to OP3. The function
+/// respects existing hydrogens when `remove_existing_h` is disabled.
+///
+/// # Arguments
+///
+/// * `residue` - Nucleic acid residue with a 5'-terminal phosphate.
+/// * `config` - Hydrogenation configuration containing pH settings.
+///
+/// # Returns
+///
+/// `Ok(())` when hydrogens are appropriately placed or removed.
+///
+/// # Errors
+///
+/// Returns [`Error::IncompleteResidueForHydro`] if phosphate atoms are missing.
+fn construct_5_prime_phosphate_hydrogens(
+    residue: &mut Residue,
+    config: &HydroConfig,
+) -> Result<(), Error> {
+    let ph = config.target_ph.unwrap_or(7.4);
+
+    if ph >= PHOSPHATE_PKA2 {
+        residue.remove_atom("HOP3");
+        residue.remove_atom("HOP2");
+        return Ok(());
+    }
+
+    if residue.has_atom("HOP3") {
+        return Ok(());
+    }
+
+    let op3 = residue
+        .atom("OP3")
+        .ok_or_else(|| Error::incomplete_for_hydro(&residue.name, residue.id, "OP3"))?
+        .pos;
+    let p = residue
+        .atom("P")
+        .ok_or_else(|| Error::incomplete_for_hydro(&residue.name, residue.id, "P"))?
+        .pos;
+
+    let direction = (op3 - p).normalize();
+    let h_pos = op3 + direction * 0.96;
+
+    residue.add_atom(Atom::new("HOP3", Element::H, h_pos));
+    Ok(())
+}
+
 /// Computes the optimal rigid transform mapping template anchor points to residue atoms.
 ///
 /// Uses Kabsch alignment with safeguards for one- and two-point configurations.
@@ -764,6 +818,64 @@ mod tests {
         residue
     }
 
+    fn five_prime_residue_with_phosphate(id: i32) -> Residue {
+        let template = db::get_template("DA").unwrap();
+        let mut residue = Residue::new(
+            id,
+            None,
+            "DA",
+            Some(StandardResidue::DA),
+            ResidueCategory::Standard,
+        );
+        residue.position = ResiduePosition::FivePrime;
+        for (atom_name, element, pos) in template.heavy_atoms() {
+            residue.add_atom(Atom::new(atom_name, element, pos));
+        }
+        let p_pos = residue.atom("P").unwrap().pos;
+        let op1_pos = residue.atom("OP1").unwrap().pos;
+        let op2_pos = residue.atom("OP2").unwrap().pos;
+        let o5_pos = residue.atom("O5'").unwrap().pos;
+        let centroid = (op1_pos.coords + op2_pos.coords + o5_pos.coords) / 3.0;
+        let direction = (p_pos.coords - centroid).normalize();
+        let op3_pos = p_pos + direction * 1.48;
+        residue.add_atom(Atom::new("OP3", Element::O, op3_pos));
+        residue
+    }
+
+    fn five_prime_residue_without_phosphate(id: i32) -> Residue {
+        let template = db::get_template("DA").unwrap();
+        let mut residue = Residue::new(
+            id,
+            None,
+            "DA",
+            Some(StandardResidue::DA),
+            ResidueCategory::Standard,
+        );
+        residue.position = ResiduePosition::FivePrime;
+        for (atom_name, element, pos) in template.heavy_atoms() {
+            if !matches!(atom_name, "P" | "OP1" | "OP2") {
+                residue.add_atom(Atom::new(atom_name, element, pos));
+            }
+        }
+        residue
+    }
+
+    fn three_prime_residue(id: i32) -> Residue {
+        let template = db::get_template("DA").unwrap();
+        let mut residue = Residue::new(
+            id,
+            None,
+            "DA",
+            Some(StandardResidue::DA),
+            ResidueCategory::Standard,
+        );
+        residue.position = ResiduePosition::ThreePrime;
+        for (atom_name, element, pos) in template.heavy_atoms() {
+            residue.add_atom(Atom::new(atom_name, element, pos));
+        }
+        residue
+    }
+
     #[test]
     fn titratable_templates_exist_in_database() {
         let expected = [
@@ -973,4 +1085,66 @@ mod tests {
         assert!(!res1.has_atom("HG"));
         assert!(!res2.has_atom("HG"));
     }
+
+    #[test]
+    fn five_prime_phosphate_deprotonated_at_physiological_ph() {
+        let residue = five_prime_residue_with_phosphate(60);
+        let mut structure = structure_with_residue(residue);
+
+        add_hydrogens(&mut structure, &HydroConfig::default()).expect("hydrogenation succeeds");
+
+        let residue = structure.find_residue("A", 60, None).unwrap();
+        assert!(residue.has_atom("OP3"), "OP3 should remain");
+        assert!(
+            !residue.has_atom("HOP3"),
+            "HOP3 should not exist at neutral pH"
+        );
+        assert!(
+            !residue.has_atom("HOP2"),
+            "HOP2 should not exist at neutral pH"
+        );
+    }
+
+    #[test]
+    fn five_prime_phosphate_protonated_below_pka() {
+        let residue = five_prime_residue_with_phosphate(61);
+        let mut structure = structure_with_residue(residue);
+        let mut config = HydroConfig::default();
+        config.target_ph = Some(5.5);
+
+        add_hydrogens(&mut structure, &config).expect("hydrogenation succeeds");
+
+        let residue = structure.find_residue("A", 61, None).unwrap();
+        assert!(residue.has_atom("OP3"), "OP3 should remain");
+        assert!(residue.has_atom("HOP3"), "HOP3 should be added below pKa");
+    }
+
+    #[test]
+    fn five_prime_without_phosphate_gets_ho5() {
+        let residue = five_prime_residue_without_phosphate(62);
+        let mut structure = structure_with_residue(residue);
+
+        add_hydrogens(&mut structure, &HydroConfig::default()).expect("hydrogenation succeeds");
+
+        let residue = structure.find_residue("A", 62, None).unwrap();
+        assert!(
+            residue.has_atom("HO5'"),
+            "HO5' should be added for 5'-OH terminus"
+        );
+        assert!(!residue.has_atom("P"), "phosphorus should not exist");
+    }
+
+    #[test]
+    fn three_prime_nucleic_gets_ho3() {
+        let residue = three_prime_residue(70);
+        let mut structure = structure_with_residue(residue);
+
+        add_hydrogens(&mut structure, &HydroConfig::default()).expect("hydrogenation succeeds");
+
+        let residue = structure.find_residue("A", 70, None).unwrap();
+        assert!(
+            residue.has_atom("HO3'"),
+            "HO3' should be added for 3' terminal"
+        );
+    }
 }