From a21b8e97a728f5e413f9b96e131e6d9cec9e6ee0 Mon Sep 17 00:00:00 2001 From: Joel Wallace Date: Sat, 15 Nov 2025 21:27:36 +0000 Subject: [PATCH] Implemented obtaining phi, psi, chi angles in main --- src/backbone/__init__.py | 0 src/backbone/ramachandran.py | 75 +++++++++++++++++++ src/legacy/backbone.py | 15 ++-- src/legacy/ramachandran.py | 57 ++++++-------- src/main.py | 15 ++++ src/rotamers/__init__.py | 0 src/rotamers/dunbrack.py | 2 +- src/rotamers/{rotamers.py => rotamers_lib.py} | 0 8 files changed, 121 insertions(+), 43 deletions(-) create mode 100644 src/backbone/__init__.py create mode 100644 src/backbone/ramachandran.py create mode 100644 src/main.py create mode 100644 src/rotamers/__init__.py rename src/rotamers/{rotamers.py => rotamers_lib.py} (100%) diff --git a/src/backbone/__init__.py b/src/backbone/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/src/backbone/ramachandran.py b/src/backbone/ramachandran.py new file mode 100644 index 0000000..fb3958a --- /dev/null +++ b/src/backbone/ramachandran.py @@ -0,0 +1,75 @@ +import numpy as np + +class RamachandranSampler: + def __init__(self): + # the weights tend towards the beta strand region to avoid collapse + + # generalised for 18 amino acids + self.general_dist = [ + # beta strand (top left) + {'phi_m': -120, 'phi_s': 20, 'psi_m': 140, 'psi_s': 20, 'w': 0.70}, + # right handed helix (bottom left) + {'phi_m': -65, 'phi_s': 15, 'psi_m': -40, 'psi_s': 15, 'w': 0.25}, + # left handed helix (top right) + {'phi_m': 60, 'phi_s': 15, 'psi_m': 40, 'psi_s': 15, 'w': 0.05} + ] + + # glycine more flexible + self.glycine_dist = [ + # beta strand (top left) + {'phi_m': -100, 'phi_s': 30, 'psi_m': 140, 'psi_s': 30, 'w': 0.3}, + # right handed helix (bottom left) + {'phi_m': -60, 'phi_s': 30, 'psi_m': -30, 'psi_s': 30, 'w': 0.2}, + # left handed helix (top right) + {'phi_m': 60, 'phi_s': 30, 'psi_m': 30, 'psi_s': 30, 'w': 0.2}, + # bottom right for glycine + {'phi_m': 100, 'phi_s': 30, 'psi_m': -140,'psi_s': 30, 'w': 0.3} + ] + + # proline no left handed helix + self.proline_dist = [ + # beta strand (top left) + {'phi_m': -63, 'phi_s': 5, 'psi_m': 150, 'psi_s': 20, 'w': 0.8}, + # right handed helix (bottom left) + {'phi_m': -63, 'phi_s': 5, 'psi_m': -35, 'psi_s': 20, 'w': 0.2} + ] + + # returns the correct distribution + def _get_distribution(self, res_name): + if res_name == 'GLY': + return self.glycine_dist + elif res_name == 'PRO': + return self.proline_dist + else: + return self.general_dist + + # returns (phi, psi) for the residue (units degrees) + def sample(self, res_name): + # the possible regions to sample + dist_options = self._get_distribution(res_name) + + weights = [d['w'] for d in dist_options] + # normalise the weights to ensure sum is exactly 1.0 (no float errors) + weights = np.array(weights) / np.sum(weights) + # choose a region based on the weights + choice_idx = np.random.choice(len(dist_options), p=weights) + selected = dist_options[choice_idx] + + # do a gaussian sample to get some variation + phi = np.random.normal(selected['phi_m'], selected['phi_s']) + psi = np.random.normal(selected['psi_m'], selected['psi_s']) + + # wrap angles to +-180 degrees + phi = (phi + 180) % 360 - 180 + psi = (psi + 180) % 360 - 180 + + return phi, psi + +# example +if __name__ == "__main__": + sampler = RamachandranSampler() + + print("Sampling 5 residues...") + for res in ['ALA', 'GLY', 'PRO', 'TRP', 'VAL']: + phi, psi = sampler.sample(res) + print(f"{res}: Phi={phi:.1f}, Psi={psi:.1f}") diff --git a/src/legacy/backbone.py b/src/legacy/backbone.py index f81167e..22c1102 100644 --- a/src/legacy/backbone.py +++ b/src/legacy/backbone.py @@ -4,23 +4,20 @@ from mpl_toolkits.mplot3d import Axes3D from nerf import place_atom_nerf from ramachandran import RamachandranSampler -# --- Constants (Idealized Geometry in Angstroms & Degrees) --- +# Ideal geometry. Angstroms for lengths, degrees for angles. GEO = { 'N_CA_len': 1.46, 'CA_C_len': 1.51, - 'C_N_len': 1.33, # Peptide bond length - - # --- MISSING KEYS ADDED HERE --- - 'N_H_len': 1.01, # Backbone Amide H bond length - 'C_N_H_angle': 119.0, # Angle for placing the H - # ------------------------------- + 'C_N_len': 1.33, + 'N_H_len': 1.01, - # Bond Angles (Standard idealized values) 'N_CA_C_angle': 111.0, 'CA_C_N_angle': 116.0, 'C_N_CA_angle': 122.0, + 'C_N_H_angle': 119.0, } +# Legacy function checking if two atoms (implementation as CAs) are within threshold. def check_clashes(new_coord, existing_coords, threshold=3.0): if len(existing_coords) == 0: return False diff = existing_coords - new_coord @@ -165,4 +162,4 @@ if __name__ == "__main__": backbone = generate_backbone_chain(test_seq, place_atom_nerf, sampler) # 4. Visualize - plot_backbone(backbone) \ No newline at end of file + plot_backbone(backbone) diff --git a/src/legacy/ramachandran.py b/src/legacy/ramachandran.py index 43087b5..58655c3 100644 --- a/src/legacy/ramachandran.py +++ b/src/legacy/ramachandran.py @@ -2,45 +2,38 @@ import numpy as np class RamachandranSampler: def __init__(self): - # Format: 'Region_Name': { - # 'phi_mean': x, 'phi_sigma': x, - # 'psi_mean': x, 'psi_sigma': x, - # 'weight': probability - # } - # 1. GENERAL (18 Amino Acids) - # Tuned for "Random Coil" (High Beta/PPII, Low Alpha) + # Generalised for 18 amino acids. Weighted for beta strand region self.general_dist = [ - # Beta / PPII Region (Favored in unfolded states) - {'phi_m': -120, 'phi_s': 20, 'psi_m': 140, 'psi_s': 20, 'w': 0.60}, - # Alpha-Right Region (Less common in unfolded, but present) - {'phi_m': -60, 'phi_s': 15, 'psi_m': -40, 'psi_s': 15, 'w': 0.30}, - # Alpha-Left / Bridge (Rare) - {'phi_m': 60, 'phi_s': 15, 'psi_m': 40, 'psi_s': 15, 'w': 0.10} + # beta strand (top left) + {'phi_m': -120, 'phi_s': 20, 'psi_m': 140, 'psi_s': 20, 'w': 0.70}, + # right handed helix (bottom left) + {'phi_m': -65, 'phi_s': 15, 'psi_m': -40, 'psi_s': 15, 'w': 0.25}, + # left handed helix (top right) + {'phi_m': 60, 'phi_s': 15, 'psi_m': 40, 'psi_s': 15, 'w': 0.05} ] - # 2. GLYCINE (Flexible) - # Glycine can visit valid regions in all 4 quadrants + # Glycine more flexible self.glycine_dist = [ - # Top Left (Beta/PPII) + # beta strand (top left) {'phi_m': -100, 'phi_s': 30, 'psi_m': 140, 'psi_s': 30, 'w': 0.3}, - # Bottom Left (Alpha-R) + # right handed helix (bottom left) {'phi_m': -60, 'phi_s': 30, 'psi_m': -30, 'psi_s': 30, 'w': 0.2}, - # Top Right (Left-handed Helix - Unique to Gly) + # left handed helix (top right) {'phi_m': 60, 'phi_s': 30, 'psi_m': 30, 'psi_s': 30, 'w': 0.2}, - # Bottom Right (Inverted Beta) + # bottom right for glycine {'phi_m': 100, 'phi_s': 30, 'psi_m': -140,'psi_s': 30, 'w': 0.3} ] - # 3. PROLINE (Rigid) - # Phi is strictly locked around -63 degrees + # Proline no left handed helix self.proline_dist = [ - # Beta/PPII (Dominant) + # beta strand (top left) {'phi_m': -63, 'phi_s': 5, 'psi_m': 150, 'psi_s': 20, 'w': 0.8}, - # Alpha (Minor) + # right handed helix (bottom left) {'phi_m': -63, 'phi_s': 5, 'psi_m': -35, 'psi_s': 20, 'w': 0.2} ] + # returns the correct distribution def _get_distribution(self, res_name): if res_name == 'GLY': return self.glycine_dist @@ -49,31 +42,29 @@ class RamachandranSampler: else: return self.general_dist - def sample(self, res_name): - """ - Returns a tuple (phi, psi) in degrees for the given residue. - """ + # returns (phi, psi) for the residue (units degrees) + def sample(self, res_name): + # the possible regions to sample dist_options = self._get_distribution(res_name) - # 1. Pick a region based on weights weights = [d['w'] for d in dist_options] - # Normalize weights to ensure sum is 1.0 (to avoid float errors) + # normalise the weights to ensure sum is exactly 1.0 (no float errors) weights = np.array(weights) / np.sum(weights) - + # choose a region based on the weights choice_idx = np.random.choice(len(dist_options), p=weights) selected = dist_options[choice_idx] - # 2. Sample Gaussian for that region + # do a gaussian sample to get some variation phi = np.random.normal(selected['phi_m'], selected['phi_s']) psi = np.random.normal(selected['psi_m'], selected['psi_s']) - # 3. Wrap angles to [-180, 180] range (optional but good practice) + # wrap angles to +-180 degrees phi = (phi + 180) % 360 - 180 psi = (psi + 180) % 360 - 180 return phi, psi -# --- Usage Example --- +# example if __name__ == "__main__": sampler = RamachandranSampler() diff --git a/src/main.py b/src/main.py new file mode 100644 index 0000000..209c53c --- /dev/null +++ b/src/main.py @@ -0,0 +1,15 @@ +from backbone.ramachandran import RamachandranSampler +from rotamers.dunbrack import DunbrackRotamerLibrary + +rs = RamachandranSampler() +rl = DunbrackRotamerLibrary() + +res_name = "PHE" + +print(f"Backbone torsion angles for {res_name}") +phi, psi = rs.sample(res_name) +print(f"phi: {phi}, psi: {psi}") + +print(f"Sidechain rotamers for {res_name}") +for rotamer in rl.rotamer_params(res_name, 100, 100): + print(rotamer.p, rotamer.chis) diff --git a/src/rotamers/__init__.py b/src/rotamers/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/src/rotamers/dunbrack.py b/src/rotamers/dunbrack.py index 2953da9..b653343 100644 --- a/src/rotamers/dunbrack.py +++ b/src/rotamers/dunbrack.py @@ -25,7 +25,7 @@ _dependent_cache = {} _independent_cache = {} -from rotamers import RotamerLibrary, RotamerParams, UnsupportedResTypeError, NoResidueRotamersError +from rotamers.rotamers_lib import RotamerLibrary, RotamerParams, UnsupportedResTypeError, NoResidueRotamersError class DunbrackRotamerLibrary(RotamerLibrary): diff --git a/src/rotamers/rotamers.py b/src/rotamers/rotamers_lib.py similarity index 100% rename from src/rotamers/rotamers.py rename to src/rotamers/rotamers_lib.py