From 13bcd4aec78cc36b60175eef486a2013a389f1cd Mon Sep 17 00:00:00 2001 From: Yi Zhou <17245097+y1zhou@users.noreply.github.com> Date: Wed, 8 May 2024 23:30:07 +0800 Subject: [PATCH] feat(rings): identify aromatic interactions (#6) * feat(ring): add fn to detect ring center and normal with SVD * feat(ring): detect rings during struct initialization * fix(ring): the U matrix should be used for finding the normal vector Our input matrix is 3xN instead of Nx3 * refactor(ring): extract center and normal into struct * feat(ring): find cation-pi interactions * feat(ring): let main fn find cation-pi interactions * feat(ring): find cation-pi and pi-pi interactions * fix(ring): wrong vectors in both angle calculations * feat(ring): use pi-pi nomenclature by Chakrabarti and Bhattacharyya * docs(readme): finish aromatic interactions --- Cargo.toml | 1 + README.md | 6 +- src/interactions/aromatic.rs | 87 +++++++++++ src/interactions/complex.rs | 284 ++++++++++++++++++++++++++--------- src/interactions/ionic.rs | 2 +- src/interactions/mod.rs | 2 + src/interactions/structs.rs | 12 +- src/main.rs | 12 +- src/residues.rs | 72 +++++++++ 9 files changed, 397 insertions(+), 81 deletions(-) create mode 100644 src/interactions/aromatic.rs diff --git a/Cargo.toml b/Cargo.toml index fc48116..a6db3b3 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -8,6 +8,7 @@ edition = "2021" [dependencies] clap = { version = "4.5.4", features = ["derive"] } +nalgebra = "0.32.5" pdbtbx = { version = "0.11.0", features = ["rayon", "rstar"] } rayon = "1.10.0" rstar = "0.12.0" diff --git a/README.md b/README.md index 7802926..7c19819 100644 --- a/README.md +++ b/README.md @@ -12,12 +12,12 @@ This is a port of the [Arpeggio](https://github.com/PDBeurope/arpeggio/) library - [x] Steric clashes - [x] VdW interactions - [x] Hydrophobic interactions - - [ ] Aromatic interactions - - [ ] Cation-pi interactions + - [x] Aromatic interactions + - [x] Cation-pi interactions - [x] Ionic interactions - [x] Hydrogen bonds - [x] Weak hydrogen bonds - - [ ] Disulfide bonds + - [x] Disulfide bonds - [ ] Covalent bonds - [ ] Output results in various formats (e.g., JSON, CSV) - [ ] Bundle into a `PyO3` extension module diff --git a/src/interactions/aromatic.rs b/src/interactions/aromatic.rs new file mode 100644 index 0000000..ff835ec --- /dev/null +++ b/src/interactions/aromatic.rs @@ -0,0 +1,87 @@ +use super::ionic::is_pos_ionizable; +use super::structs::Interaction; +use crate::residues::Ring; + +use nalgebra as na; +use pdbtbx::*; + +const CATION_PI_ANGLE_THRESHOLD: f64 = 30.0; +const CATION_PI_DIST_THRESHOLD: f64 = 4.5; +const PI_PI_DIST_THRESHOLD: f64 = 6.0; + +/// Identify cation-pi interactions. +pub fn find_cation_pi(ring: &Ring, entity: &AtomConformerResidueChainModel) -> Option { + if is_pos_ionizable(entity.residue().name().unwrap(), entity.atom().name()) { + let atom_coord = entity.atom().pos(); + let atom_point = na::Vector3::new(atom_coord.0, atom_coord.1, atom_coord.2); + let dist = point_ring_dist(ring, &atom_coord); + let theta = point_ring_angle(ring, &atom_point); + + if (theta <= CATION_PI_ANGLE_THRESHOLD) & (dist <= CATION_PI_DIST_THRESHOLD) { + Some(Interaction::CationPi) + } else { + None + } + } else { + None + } +} + +/// Identify pi-pi interactions using the classification by [Chakrabarti and Bhattacharyya (2007)](https://doi.org/10.1016/j.pbiomolbio.2007.03.016). +pub fn find_pi_pi(ring1: &Ring, ring2: &Ring) -> Option { + let angle_vec = ring1.center - ring2.center; + let dist = (angle_vec).norm(); + if dist <= PI_PI_DIST_THRESHOLD { + let theta = point_ring_angle(ring1, &angle_vec); + let dihedral = ring_ring_angle(ring1, ring2); + + match dihedral { + d if d <= 30.0 => match theta { + t if t <= 30.0 => Some(Interaction::PiSandwichStacking), // ff + t if t <= 60.0 => Some(Interaction::PiDisplacedStacking), // of + t if t <= 90.0 => Some(Interaction::PiParallelInPlaneStacking), // ee + _ => None, + }, + // ft, ot, and et + d if d <= 60.0 => Some(Interaction::PiTiltedStacking), + d if d <= 90.0 => match theta { + t if (30.0..60.0).contains(&t) => Some(Interaction::PiLStacking), // oe + _ => Some(Interaction::PiTStacking), // fe and ef + }, + _ => None, + } + } else { + None + } +} + +/// Calculate the distance from the point to the ring center. +pub fn point_ring_dist(ring: &Ring, point: &(f64, f64, f64)) -> f64 { + let atom_point = na::Vector3::new(point.0, point.1, point.2); + (atom_point - ring.center).norm() +} + +/// Calculate the angle between the ring normal and the vector pointing from +/// the ring center to the point. +fn point_ring_angle(ring: &Ring, point: &na::Vector3) -> f64 { + let v = point - ring.center; + let mut rad = (ring.normal.dot(&v) / (ring.normal.norm() * v.norm())).acos(); + + // Convert to degrees + if rad > std::f64::consts::FRAC_PI_2 { + rad = std::f64::consts::PI - rad; + } + rad.to_degrees() +} + +/// Calculate the angle between two rings. +fn ring_ring_angle(ring1: &Ring, ring2: &Ring) -> f64 { + let mut rad = + (ring1.normal.dot(&ring2.normal) / (ring1.normal.norm() * ring2.normal.norm())).acos(); + + // Convert to degrees + if rad > std::f64::consts::FRAC_PI_2 { + rad = std::f64::consts::PI - rad; + } + rad.to_degrees() +} diff --git a/src/interactions/complex.rs b/src/interactions/complex.rs index 4056ea9..05ac275 100644 --- a/src/interactions/complex.rs +++ b/src/interactions/complex.rs @@ -1,8 +1,12 @@ use super::{ - find_hydrogen_bond, find_hydrophobic_contact, find_ionic_bond, find_vdw_contact, - find_weak_hydrogen_bond, Interaction, ResultEntry, + find_cation_pi, find_hydrogen_bond, find_hydrophobic_contact, find_ionic_bond, find_pi_pi, + find_vdw_contact, find_weak_hydrogen_bond, point_ring_dist, InteractingEntity, Interaction, + ResultEntry, +}; +use crate::{ + residues::{ResidueExt, Ring}, + utils::{hierarchy_to_entity, parse_groups}, }; -use crate::utils::{hierarchy_to_entity, parse_groups}; use pdbtbx::*; use rayon::prelude::*; @@ -23,15 +27,22 @@ pub struct InteractionComplex { /// Maps residue names to unique indices res2idx: HashMap>, + /// Maps ring residues to ring centers and normals + rings: HashMap<(String, isize, String, String), Ring>, } impl InteractionComplex { pub fn new(model: PDB, groups: &str, vdw_comp_factor: f64, interacting_threshold: f64) -> Self { + // Parse all chains and input chain groups let all_chains: HashSet = model.par_chains().map(|c| c.id().to_string()).collect(); let (ligand, receptor) = parse_groups(&all_chains, groups); + // Build a mapping of residue names to indices let res2idx = build_residue_index(&model); + // Build a mapping of ring residue names to ring centers and normals + let rings = build_ring_positions(&model); + Self { model, ligand, @@ -39,51 +50,62 @@ impl InteractionComplex { vdw_comp_factor, interacting_threshold, res2idx, + rings, } } - fn is_neighboring_res_pair( + fn is_neighboring_hierarchy( &self, x: &AtomConformerResidueChainModel, y: &AtomConformerResidueChainModel, ) -> bool { - if x.chain().id() != y.chain().id() { - return false; - } let (x_resi, x_insertion) = x.residue().id(); - let x_altloc = match x_insertion { - Some(insertion) => insertion.to_string(), - None => "".to_string(), - }; - let x_idx = self.res2idx[x.chain().id()][&(x_resi, x_altloc)]; + let x_altloc = x_insertion.unwrap_or(""); let (y_resi, y_insertion) = y.residue().id(); - let y_altloc = match y_insertion { - Some(insertion) => insertion.to_string(), - None => "".to_string(), - }; - let y_idx = self.res2idx[y.chain().id()][&(y_resi, y_altloc)]; + let y_altloc = y_insertion.unwrap_or(""); + self.is_neighboring_res_pair( + x.chain().id(), + x_resi, + x_altloc, + y.chain().id(), + y_resi, + y_altloc, + ) + } + + fn is_neighboring_res_pair( + &self, + x_chain: &str, + x_resi: isize, + x_altloc: &str, + y_chain: &str, + y_resi: isize, + y_altloc: &str, + ) -> bool { + if x_chain != y_chain { + return false; + } + let x_idx = self.res2idx[&x_chain.to_string()][&(x_resi, x_altloc.to_string())]; + let y_idx = self.res2idx[&y_chain.to_string()][&(y_resi, y_altloc.to_string())]; match x_idx { - 0 => y_idx == x_idx + 1, - _ => (y_idx == x_idx - 1) | (y_idx == x_idx + 1), + 0 => (y_idx == x_idx) | (y_idx == x_idx + 1), + _ => (y_idx == x_idx - 1) | (y_idx == x_idx) | (y_idx == x_idx + 1), } } } pub trait Interactions { - /// Get all protein-protein interactions between the ligand and receptor - fn get_ppi(&self); - /// Get all atomic interactions between the ligand and receptor + /// Get all atomic interactions between the ligand and receptor. fn get_atomic_contacts(&self) -> Vec; - /// Get all ring interactions between the ligand and receptor - fn get_ring_contacts(&self); + + /// Get all ring-atom interactions between the ligand and receptor. + fn get_ring_atom_contacts(&self) -> Vec; + /// Get all ring-ring interactions between the ligand and receptor. + fn get_ring_ring_contacts(&self) -> Vec; } impl Interactions for InteractionComplex { - fn get_ppi(&self) { - todo!() - } - fn get_atomic_contacts(&self) -> Vec { let tree = self.model.create_hierarchy_rtree(); let max_radius_squared = self.interacting_threshold * self.interacting_threshold; @@ -91,7 +113,7 @@ impl Interactions for InteractionComplex { // Find all atoms within the radius of the ligand atoms let ligand_neighbors: Vec<( AtomConformerResidueChainModel, - AtomConformerResidueChainModel, + &AtomConformerResidueChainModel, )> = self .model .atoms_with_hierarchy() @@ -99,11 +121,7 @@ impl Interactions for InteractionComplex { self.ligand.contains(x.chain().id()) & (x.atom().element().unwrap() != &Element::H) }) .flat_map(|x| { - let neighbor_entities: Vec<( - AtomConformerResidueChainModel, - AtomConformerResidueChainModel, - )> = tree - .locate_within_distance(x.atom().pos(), max_radius_squared) + tree.locate_within_distance(x.atom().pos(), max_radius_squared) .filter(|y| { self.receptor.contains(y.chain().id()) & (y.atom().element().unwrap() != &Element::H) @@ -111,35 +129,28 @@ impl Interactions for InteractionComplex { & !((x.chain().id() == y.chain().id()) & (x.residue().serial_number() > y.residue().serial_number())) // Skip neighboring residues - & !self.is_neighboring_res_pair(&x, y) - // Skip atoms in the same residue - & !((x.chain().id() == y.chain().id()) - & (x.residue().serial_number() == y.residue().serial_number()) - & (x.residue().insertion_code() == y.residue().insertion_code())) + & !self.is_neighboring_hierarchy(&x, y) }) - .map(|y| (x.clone(), y.clone())) - .collect(); - - match neighbor_entities.len() { - 0 => vec![], - _ => neighbor_entities, - } + .map(|y| (x.clone(), y)) + .collect::>() }) .collect(); ligand_neighbors .par_iter() - .filter_map(|x| { + .filter_map(|(e1, e2)| { let mut atomic_contacts: Vec = vec![]; // Clashes and VdW contacts - let vdw = - find_vdw_contact(&x.0, &x.1, self.vdw_comp_factor).map(|intxn| ResultEntry { - interaction: intxn, - ligand: hierarchy_to_entity(&x.0), - receptor: hierarchy_to_entity(&x.1), - distance: x.0.atom().distance(x.1.atom()), - }); + let vdw = find_vdw_contact(e1, e2, self.vdw_comp_factor).map(|intxn| ResultEntry { + interaction: intxn, + ligand: hierarchy_to_entity(e1), + receptor: hierarchy_to_entity(e2), + distance: e1.atom().distance(e2.atom()), + }); atomic_contacts.extend(vdw.clone()); // Skip checking for other interactions if there is a clash @@ -149,42 +160,42 @@ impl Interactions for InteractionComplex { // Hydrogen bonds and polar contacts let hbonds = - find_hydrogen_bond(&x.0, &x.1, self.vdw_comp_factor).map(|intxn| ResultEntry { + find_hydrogen_bond(e1, e2, self.vdw_comp_factor).map(|intxn| ResultEntry { interaction: intxn, - ligand: hierarchy_to_entity(&x.0), - receptor: hierarchy_to_entity(&x.1), - distance: x.0.atom().distance(x.1.atom()), + ligand: hierarchy_to_entity(e1), + receptor: hierarchy_to_entity(e2), + distance: e1.atom().distance(e2.atom()), }); atomic_contacts.extend(hbonds); // C-H...O bonds let weak_hbonds = - find_weak_hydrogen_bond(&x.0, &x.1, self.vdw_comp_factor).map(|intxn| { + find_weak_hydrogen_bond(e1, e2, self.vdw_comp_factor).map(|intxn| { ResultEntry { interaction: intxn, - ligand: hierarchy_to_entity(&x.0), - receptor: hierarchy_to_entity(&x.1), - distance: x.0.atom().distance(x.1.atom()), + ligand: hierarchy_to_entity(e1), + receptor: hierarchy_to_entity(e2), + distance: e1.atom().distance(e2.atom()), } }); atomic_contacts.extend(weak_hbonds); // Ionic bonds - let ionic_bonds = find_ionic_bond(&x.0, &x.1).map(|intxn| ResultEntry { + let ionic_bonds = find_ionic_bond(e1, e2).map(|intxn| ResultEntry { interaction: intxn, - ligand: hierarchy_to_entity(&x.0), - receptor: hierarchy_to_entity(&x.1), - distance: x.0.atom().distance(x.1.atom()), + ligand: hierarchy_to_entity(e1), + receptor: hierarchy_to_entity(e2), + distance: e1.atom().distance(e2.atom()), }); atomic_contacts.extend(ionic_bonds); // Hydrophobic contacts let hydrophobic_contacts = - find_hydrophobic_contact(&x.0, &x.1).map(|intxn| ResultEntry { + find_hydrophobic_contact(e1, e2).map(|intxn| ResultEntry { interaction: intxn, - ligand: hierarchy_to_entity(&x.0), - receptor: hierarchy_to_entity(&x.1), - distance: x.0.atom().distance(x.1.atom()), + ligand: hierarchy_to_entity(e1), + receptor: hierarchy_to_entity(e2), + distance: e1.atom().distance(e2.atom()), }); atomic_contacts.extend(hydrophobic_contacts); @@ -194,8 +205,112 @@ impl Interactions for InteractionComplex { .collect::>() } - fn get_ring_contacts(&self) { - todo!() + fn get_ring_atom_contacts(&self) -> Vec { + let tree = self.model.create_hierarchy_rtree(); + let max_radius_squared = self.interacting_threshold * self.interacting_threshold; + + // Find ring - atom contacts + let ring_atom_neighbors = self + .rings + .iter() + .flat_map(|(k, v)| { + tree + .locate_within_distance((v.center.x, v.center.y, v.center.z), max_radius_squared) + // Ring on ligand or receptor, and atom on the other side + .filter(|y| + (self.ligand.contains(k.0.as_str()) + & self.receptor.contains(y.chain().id()) + ) | + (self.ligand.contains(y.chain().id()) + & self.receptor.contains(k.0.as_str()) + ) + ) + .filter(|y| + (y.atom().element().unwrap() != &Element::H) + // Skip lower resi if both entities are on the same chain + & !((k.0 == y.chain().id()) + & (k.1 > y.residue().serial_number())) + // Skip neighboring residues + & !self.is_neighboring_res_pair(k.0.as_str(), k.1, k.2.as_str(), y.chain().id(), y.residue().serial_number(), y.residue().insertion_code().unwrap_or("")) + ).map(|y| (k, v, y)) + .collect::>() + }).collect::>(); + + // Find ring-atom interactions + ring_atom_neighbors + .par_iter() + .filter_map(|(k, ring, x)| { + let mut ring_contacts = Vec::new(); + + // Cation-pi interactions + let dist = point_ring_dist(ring, &x.atom().pos()); + let cation_pi_contacts = find_cation_pi(ring, x).map(|intxn| ResultEntry { + interaction: intxn, + ligand: InteractingEntity { + chain: k.0.clone(), + resi: k.1, + altloc: k.2.clone(), + resn: k.3.clone(), + atomn: "Ring".to_string(), + atomi: 0, + }, + receptor: hierarchy_to_entity(x), + distance: dist, + }); + ring_contacts.extend(cation_pi_contacts); + + Some(ring_contacts) + }) + .flatten() + .collect::>() + } + + fn get_ring_ring_contacts(&self) -> Vec { + // Find ring - ring contacts + let ring_ring_neighbors = self + .rings + .iter() + .filter(|(k, _)| self.ligand.contains(k.0.as_str())) + .flat_map(|(k1, ring1)| { + self.rings.iter() + .filter(|(k2, _)| + self.receptor.contains(k2.0.as_str()) + & !((k1.0 == k2.0) & (k1.1 > k2.1)) // Skip lower resi if both entities are on the same chain + & !self.is_neighboring_res_pair(k1.0.as_str(), k1.1, k1.2.as_str(), k2.0.as_str(), k2.1, k2.2.as_str()) // Skip neighboring residues + ).map(|(k2, ring2)| (k1, ring1, k2, ring2)) + .collect::>() + }).collect::>(); + + // Find ring-ring interactions + ring_ring_neighbors + .par_iter() + .filter_map(|(k1, ring1, k2, ring2)| { + let dist = (ring1.center - ring2.center).norm(); + let pi_pi_contacts = find_pi_pi(ring1, ring2).map(|intxn| ResultEntry { + interaction: intxn, + ligand: InteractingEntity { + chain: k1.0.clone(), + resi: k1.1, + altloc: k1.2.clone(), + resn: k1.3.clone(), + atomn: "Ring".to_string(), + atomi: 0, + }, + receptor: InteractingEntity { + chain: k2.0.clone(), + resi: k2.1, + altloc: k2.2.clone(), + resn: k2.3.clone(), + atomn: "Ring".to_string(), + atomi: 0, + }, + distance: dist, + }); + + Some(pi_pi_contacts) + }) + .flatten() + .collect::>() } } @@ -225,3 +340,28 @@ fn build_residue_index(model: &PDB) -> HashMap HashMap<(String, isize, String, String), Ring> { + let ring_res = HashSet::from(["HIS", "PHE", "TYR", "TRP"]); + model + .atoms_with_hierarchy() + .filter(|x| ring_res.contains(x.residue().name().unwrap())) + .map(|x| { + let (resi, insertion) = x.residue().id(); + let altloc = match insertion { + Some(insertion) => insertion.to_string(), + None => "".to_string(), + }; + + ( + ( + x.chain().id().to_string(), + resi, + altloc, + x.residue().name().unwrap().to_string(), + ), + x.residue().ring_center_and_normal().unwrap(), + ) + }) + .collect::>() +} diff --git a/src/interactions/ionic.rs b/src/interactions/ionic.rs index 665b5e1..be0b3c3 100644 --- a/src/interactions/ionic.rs +++ b/src/interactions/ionic.rs @@ -47,7 +47,7 @@ fn is_ionic_pair<'a>( } /// Check if the entity contains ionizable groups that are positively charged at pH 7.0. -fn is_pos_ionizable(res_name: &str, atom_name: &str) -> bool { +pub fn is_pos_ionizable(res_name: &str, atom_name: &str) -> bool { matches!( (res_name, atom_name), ("ARG", "NE") diff --git a/src/interactions/mod.rs b/src/interactions/mod.rs index 9ba893d..0a7ab39 100644 --- a/src/interactions/mod.rs +++ b/src/interactions/mod.rs @@ -1,3 +1,4 @@ +pub mod aromatic; pub mod complex; pub mod hbond; pub mod hydrophobic; @@ -6,6 +7,7 @@ pub mod structs; pub mod vdw; // Re-exports +pub use aromatic::{find_cation_pi, find_pi_pi, point_ring_dist}; pub use complex::*; pub use hbond::{find_hydrogen_bond, find_weak_hydrogen_bond}; pub use hydrophobic::find_hydrophobic_contact; diff --git a/src/interactions/structs.rs b/src/interactions/structs.rs index 0df73b8..847c826 100644 --- a/src/interactions/structs.rs +++ b/src/interactions/structs.rs @@ -27,12 +27,18 @@ pub enum Interaction { WeakPolarContact, // Aromatic - /// Pi-pi staggered stacking, parallel displaced + /// Pi-pi staggered stacking, parallel displaced (`of`) PiDisplacedStacking, - /// Pi-pi perpendicular T-shaped stacking + /// Pi-pi perpendicular T-shaped stacking (`fe` and `ef`) PiTStacking, - /// Pi-pi direct stacking, repulsive + /// Pi-pi direct stacking, repulsive (`ff`) PiSandwichStacking, + /// Pi-pi parallel in-plane interactions (`ee`) + PiParallelInPlaneStacking, + /// Pi-pi tilted interactions (`ft`, `ot`, and `et`) + PiTiltedStacking, + /// Pi-pi cogwheel or L-shaped interactions (`oe`) + PiLStacking, /// Cation-pi interaction CationPi, diff --git a/src/main.rs b/src/main.rs index c4e1fab..b2da115 100644 --- a/src/main.rs +++ b/src/main.rs @@ -1,4 +1,5 @@ #![warn(missing_docs)] +#![allow(clippy::type_complexity)] #![doc = include_str!("../README.md")] mod chains; @@ -7,7 +8,7 @@ mod residues; mod utils; use crate::chains::ChainExt; -use crate::interactions::{Interaction, InteractionComplex, Interactions}; +use crate::interactions::{Interaction, InteractionComplex, Interactions, ResultEntry}; use crate::utils::load_model; use clap::Parser; @@ -88,6 +89,7 @@ fn main() { receptor = i_complex.receptor ); + // Find interactions let atomic_contacts = i_complex.get_atomic_contacts(); info!("Found {} atom-atom contacts", atomic_contacts.len()); atomic_contacts.iter().for_each(|h| { @@ -95,5 +97,11 @@ fn main() { Interaction::StericClash => warn!("{h}"), _ => debug!("{h}"), }; - }) + }); + + let mut ring_contacts: Vec = Vec::new(); + ring_contacts.extend(i_complex.get_ring_atom_contacts()); + ring_contacts.extend(i_complex.get_ring_ring_contacts()); + info!("Found {} ring contacts", ring_contacts.len()); + ring_contacts.iter().for_each(|h| debug!("{h}")); } diff --git a/src/residues.rs b/src/residues.rs index 1b79aef..3f75542 100644 --- a/src/residues.rs +++ b/src/residues.rs @@ -1,7 +1,20 @@ +use nalgebra as na; use pdbtbx::*; +use rayon::prelude::*; + +#[derive(Clone, Copy, Debug)] +pub struct Ring { + pub center: na::Vector3, + pub normal: na::Vector3, +} pub trait ResidueExt { + /// The residue one-letter code, or `None` if it's not an amino acid. fn resn(&self) -> Option<&str>; + /// Return the atoms in the aromatic ring of the residue. + fn ring_atoms(&self) -> Vec<&Atom>; + /// + fn ring_center_and_normal(&self) -> Option; } impl ResidueExt for Residue { @@ -36,4 +49,63 @@ impl ResidueExt for Residue { _ => Some(aa_code), } } + + fn ring_atoms(&self) -> Vec<&Atom> { + let res_name = self + .name() + .unwrap_or("Some conformers have different names!"); + + match res_name { + "HIS" => self + .par_atoms() + .filter(|atom| matches!(atom.name(), "CG" | "ND1" | "CE1" | "NE2" | "CD2")) + .collect(), + "PHE" | "TYR" => self + .par_atoms() + .filter(|atom| matches!(atom.name(), "CG" | "CD1" | "CD2" | "CE1" | "CE2" | "CZ")) + .collect(), + "TRP" => self + .par_atoms() + .filter(|atom| { + matches!( + atom.name(), + "CG" | "CD1" | "CD2" | "NE1" | "CE2" | "CE3" | "CZ2" | "CZ3" | "CH2" + ) + }) + .collect(), + _ => vec![], + } + } + + fn ring_center_and_normal(&self) -> Option { + let ring_atoms = self.ring_atoms(); + if ring_atoms.is_empty() { + return None; + } + + // Construct 3*N matrix of ring atom coordinates + let mut atom_coords: na::Matrix3xX = na::Matrix3xX::::from_iterator( + ring_atoms.len(), + ring_atoms.iter().flat_map(|atom| { + let coord = atom.pos(); + [coord.0, coord.1, coord.2].into_iter() + }), + ); + let center = atom_coords.column_mean(); + + // Center the matrix and perform SVD + // Ref: https://en.wikipedia.org/wiki/Singular_value_decomposition#Total_least_squares_minimization + // Ref: https://math.stackexchange.com/a/2810220 + for i in 0..atom_coords.ncols() { + atom_coords.set_column(i, &(atom_coords.column(i) - center)); + } + + let svd = atom_coords.svd(true, true); + let normal = svd.u.unwrap().column(2).clone_owned(); + // Sanity check the normal is orthogonal to the plane vectors + // let dot_products = atom_coords.transpose() * normal; + // debug!("Dot products:\n{:?}", dot_products); + + Some(Ring { center, normal }) + } }