diff --git a/MANIFEST.in b/MANIFEST.in index 0a189be..fdfa9f5 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -4,3 +4,6 @@ include propka.cfg include versioneer.py include propka/_version.py + +# Python type annotation declarations, see PEP-561 +include propka/py.typed diff --git a/propka/atom.py b/propka/atom.py index 5cb6a4e..0cca0e1 100644 --- a/propka/atom.py +++ b/propka/atom.py @@ -8,6 +8,7 @@ import string from typing import cast, List, NoReturn, Optional, TYPE_CHECKING +import warnings from propka.lib import make_tidy_atom_label from . import hybrid36 @@ -42,6 +43,43 @@ class Atom: :meth:`make_input_line` and :meth:`get_input_parameters` have been removed as reading/writing PROPKA input is no longer supported. """ + group: Optional["Group"] = None + group_type: Optional[str] = None + cysteine_bridge: bool = False + residue: NoReturn = None # type: ignore[assignment] + conformation_container: Optional["ConformationContainer"] = None + molecular_container: Optional["MolecularContainer"] = None + is_protonated: bool = False + steric_num_lone_pairs_set: bool = False + terminal: Optional[str] = None + charge: float = 0.0 + charge_set: bool = False + steric_number: int = 0 + number_of_lone_pairs: int = 0 + number_of_protons_to_add: int = 0 + num_pi_elec_2_3_bonds: int = 0 + num_pi_elec_conj_2_3_bonds: int = 0 + groups_extracted: bool = False + + # PDB attributes + name: str = '' + numb: int = 0 + x: float = 0.0 + y: float = 0.0 + z: float = 0.0 + res_num: int = 0 + res_name: str = '' + chain_id: str = 'A' + type: str = '' + occ: str = '1.0' + beta: str = '0.0' + element: str = '' + icode: str = '' + + # ligand atom types + sybyl_type = '' + sybyl_assigned = False + marvin_pka = False def __init__(self, line: Optional[str] = None): """Initialize Atom object. @@ -50,53 +88,17 @@ def __init__(self, line: Optional[str] = None): line: Line from a PDB file to set properties of atom. """ self.number_of_bonded_elements: NoReturn = cast(NoReturn, {}) # FIXME unused? - self.group: Optional[Group] = None - self.group_type: Optional[str] = None - self.cysteine_bridge: bool = False self.bonded_atoms: List[Atom] = [] - self.residue = None - self.conformation_container: Optional[ConformationContainer] = None - self.molecular_container: Optional[MolecularContainer] = None - self.is_protonated = False - self.steric_num_lone_pairs_set = False - self.terminal: Optional[str] = None - self.charge = 0.0 - self.charge_set = False - self.steric_number = 0 - self.number_of_lone_pairs = 0 - self.number_of_protons_to_add = 0 - self.num_pi_elec_2_3_bonds = 0 - self.num_pi_elec_conj_2_3_bonds = 0 - self.groups_extracted = 0 self.set_properties(line) fmt = "{r.name:3s}{r.res_num:>4d}{r.chain_id:>2s}" self.residue_label = fmt.format(r=self) - # ligand atom types - self.sybyl_type = '' - self.sybyl_assigned = False - self.marvin_pka = False - def set_properties(self, line: Optional[str]): """Line from PDB file to set properties of atom. Args: line: PDB file line """ - self.name = '' - self.numb = 0 - self.x = 0.0 - self.y = 0.0 - self.z = 0.0 - self.res_num = 0 - self.res_name = '' - self.chain_id = 'A' - self.type = '' - self.occ = '1.0' - self.beta = '0.0' - self.element = '' - self.icode = '' - if line: self.name = line[12:16].strip() self.numb = int(hybrid36.decode(line[6:11])) @@ -183,9 +185,17 @@ def is_atom_within_bond_distance(self, other_atom, max_bonds, cur_bond): return True return False - def set_property(self, numb=None, name=None, res_name=None, chain_id=None, - res_num=None, x=None, y=None, z=None, occ=None, - beta=None): + def set_property(self, + numb: Optional[int] = None, + name: Optional[str] = None, + res_name: Optional[str] = None, + chain_id: Optional[str] = None, + res_num: Optional[int] = None, + x: Optional[float] = None, + y: Optional[float] = None, + z: Optional[float] = None, + occ: Optional[str] = None, + beta: Optional[str] = None): """Set properties of the atom object. Args: @@ -303,6 +313,7 @@ def make_pdb_line2(self, numb=None, name=None, res_name=None, chain_id=None, Returns: String with PDB line. """ + warnings.warn("only used by unused function") if numb is None: numb = self.numb if name is None: @@ -343,11 +354,12 @@ def __str__(self): """Return an undefined-format string version of this atom.""" return STR_FMT.format(r=self) - def set_residue(self, residue): + def set_residue(self, residue: NoReturn): """ Makes a reference to the parent residue Args: residue: the parent residue """ + raise NotImplementedError("unused") if self.residue is None: self.residue = residue diff --git a/propka/bonds.py b/propka/bonds.py index 396ec86..9bc00c0 100644 --- a/propka/bonds.py +++ b/propka/bonds.py @@ -93,6 +93,7 @@ def find_bonds_for_protein(self, protein): Args: protein: the protein to search for bonds """ + raise NotImplementedError("unused") _LOGGER.info('++++ Side chains ++++') # side chains for chain in protein.chains: @@ -132,6 +133,7 @@ def check_for_cysteine_bonds(self, cys1, cys2): cys1: one of the cysteines to check cys1: one of the cysteines to check """ + raise NotImplementedError("unused") for atom1 in cys1.atoms: if atom1.name == 'SG': for atom2 in cys2.atoms: diff --git a/propka/conformation_container.py b/propka/conformation_container.py index 0196b99..27601c4 100644 --- a/propka/conformation_container.py +++ b/propka/conformation_container.py @@ -6,7 +6,10 @@ """ import logging import functools -from typing import Iterable, List, NoReturn, Optional, TYPE_CHECKING, Set +from typing import Callable, Dict, Iterable, Iterator, List, NoReturn, Optional, TYPE_CHECKING, Set + +from propka.lib import Options +from propka.version import Version if TYPE_CHECKING: from propka.atom import Atom @@ -19,10 +22,13 @@ from propka.determinants import set_backbone_determinants, set_ion_determinants from propka.determinants import set_determinants from propka.group import Group, is_group +from propka.parameters import Parameters _LOGGER = logging.getLogger(__name__) +CallableGroupToGroups = Callable[[Group], List[Group]] + #: A large number that gets multipled with the integer obtained from applying #: :func:`ord` to the atom chain ID. Used in calculating atom keys for @@ -44,9 +50,9 @@ class ConformationContainer: """ def __init__(self, - name: str = '', - parameters=None, - molecular_container: Optional["MolecularContainer"] = None): + name: str, + parameters: Parameters, + molecular_container: "MolecularContainer"): """Initialize conformation container. Args: @@ -145,6 +151,7 @@ def find_bonded_titratable_groups(self, atom: "Atom", num_bonds: int, Returns: a set of bonded atom groups """ + assert self.parameters is not None res: Set[Group] = set() for bond_atom in atom.bonded_atoms: # skip the original atom @@ -187,7 +194,7 @@ def init_group(self, group: Group): # If --titrate_only option is set, make non-specified residues # un-titratable: - assert self.molecular_container is not None + assert self.molecular_container.options is not None titrate_only = self.molecular_container.options.titrate_only if titrate_only is not None: atom = group.atom @@ -196,7 +203,7 @@ def init_group(self, group: Group): if group.residue_type == 'CYS': group.exclude_cys_from_results = True - def calculate_pka(self, version, options): + def calculate_pka(self, version: Version, options: Options): """Calculate pKas for conformation container. Args: @@ -272,7 +279,7 @@ def coupling_effects(self): return penalised_labels @staticmethod - def share_determinants(groups): + def share_determinants(groups: Iterable[Group]): """Share sidechain, backbone, and Coloumb determinants between groups. Args: @@ -282,7 +289,7 @@ def share_determinants(groups): types = ['sidechain', 'backbone', 'coulomb'] for type_ in types: # find maximum value for each determinant - max_dets = {} + max_dets: Dict[Group, float] = {} for group in groups: for det in group.determinants[type_]: # update max dets @@ -298,7 +305,11 @@ def share_determinants(groups): for group in groups: group.set_determinant(new_determinant, type_) - def get_coupled_systems(self, groups, get_coupled_groups): + def get_coupled_systems( + self, + groups: Iterable[Group], + get_coupled_groups: CallableGroupToGroups, + ) -> Iterator[Set[Group]]: """A generator that yields covalently coupled systems. Args: @@ -310,15 +321,16 @@ def get_coupled_systems(self, groups, get_coupled_groups): groups = set(groups) while len(groups) > 0: # extract a system of coupled groups ... - system = set() + system: Set[Group] = set() self.get_a_coupled_system_of_groups( groups.pop(), system, get_coupled_groups) # ... and remove them from the list groups -= system yield system - def get_a_coupled_system_of_groups(self, new_group, coupled_groups, - get_coupled_groups): + def get_a_coupled_system_of_groups(self, new_group: Group, + coupled_groups: Set[Group], + get_coupled_groups: CallableGroupToGroups): """Set up coupled systems of groups. Args: @@ -349,7 +361,7 @@ def calculate_folding_energy(self, ph=None, reference=None): reference=reference) return ddg - def calculate_charge(self, parameters, ph=None): + def calculate_charge(self, parameters: Parameters, ph: float): """Calculate charge for folded and unfolded states. Args: @@ -367,7 +379,7 @@ def calculate_charge(self, parameters, ph=None): state='folded') return unfolded, folded - def get_backbone_groups(self): + def get_backbone_groups(self) -> List[Group]: """Get backbone groups needed for the pKa calculations. Returns: diff --git a/propka/coupled_groups.py b/propka/coupled_groups.py index 9d9e137..7ed0e87 100644 --- a/propka/coupled_groups.py +++ b/propka/coupled_groups.py @@ -6,9 +6,11 @@ """ import logging import itertools +from typing import Optional import propka.lib from propka.group import Group from propka.output import make_interaction_map +from propka.parameters import Parameters _LOGGER = logging.getLogger(__name__) @@ -16,9 +18,8 @@ class NonCovalentlyCoupledGroups: """Groups that are coupled without covalent bonding.""" - def __init__(self): - self.parameters = None - self.do_prot_stat = True + parameters: Optional[Parameters] = None + do_prot_stat = True def is_coupled_protonation_state_probability(self, group1, group2, energy_method, @@ -33,6 +34,7 @@ def is_coupled_protonation_state_probability(self, group1, group2, Returns: dictionary describing coupling """ + assert self.parameters is not None # check if the interaction energy is high enough interaction_energy = max(self.get_interaction(group1, group2), self.get_interaction(group2, group1)) @@ -105,6 +107,7 @@ def get_pka_diff_factor(self, pka1, pka2): Returns: float value of scaling factor """ + assert self.parameters is not None intrinsic_pka_diff = abs(pka1-pka2) res = 0.0 if intrinsic_pka_diff <= self.parameters.max_intrinsic_pka_diff: @@ -122,6 +125,7 @@ def get_free_energy_diff_factor(self, energy1, energy2): Returns: float value of scaling factor """ + assert self.parameters is not None free_energy_diff = abs(energy1-energy2) res = 0.0 if free_energy_diff <= self.parameters.max_free_energy_diff: @@ -136,6 +140,7 @@ def get_interaction_factor(self, interaction_energy): Returns: float value of scaling factor """ + assert self.parameters is not None res = 0.0 interaction_energy = abs(interaction_energy) if interaction_energy >= self.parameters.min_interaction_energy: @@ -260,7 +265,7 @@ def print_system(self, conformation, system): _LOGGER.info(swap_info) @staticmethod - def get_interaction(group1, group2, include_side_chain_hbs=True): + def get_interaction(group1: Group, group2: Group, include_side_chain_hbs=True): """Get interaction energy between two groups. Args: diff --git a/propka/determinant.py b/propka/determinant.py index 03f06d7..741390f 100644 --- a/propka/determinant.py +++ b/propka/determinant.py @@ -15,6 +15,11 @@ """ +from typing import TYPE_CHECKING + +if TYPE_CHECKING: + from propka.group import Group + class Determinant: """Determinant class. @@ -25,7 +30,7 @@ class Determinant: TODO - figure out what this class does. """ - def __init__(self, group, value): + def __init__(self, group: "Group", value: float): """Initialize the object. Args: @@ -36,7 +41,7 @@ def __init__(self, group, value): self.label = group.label self.value = value - def add(self, value): + def add(self, value: float): """Increment determinant value. Args: diff --git a/propka/determinants.py b/propka/determinants.py index 235ab43..320defe 100644 --- a/propka/determinants.py +++ b/propka/determinants.py @@ -14,12 +14,18 @@ """ import math +from typing import List + +import propka.calculations import propka.iterative import propka.lib import propka.vector_algebra from propka.calculations import squared_distance, get_smallest_distance from propka.energy import angle_distance_factors, hydrogen_bond_energy from propka.determinant import Determinant +from propka.group import Group +from propka.iterative import Interaction +from propka.version import Version # Cutoff for angle factor @@ -28,7 +34,7 @@ FANGLE_MIN = 0.001 -def set_determinants(propka_groups, version=None, options=None): +def set_determinants(propka_groups: List[Group], version: Version, options=None): """Add side-chain and coulomb determinants/perturbations to all residues. NOTE - backbone determinants are set separately @@ -38,7 +44,7 @@ def set_determinants(propka_groups, version=None, options=None): version: version object options: options object """ - iterative_interactions = [] + iterative_interactions: List[Interaction] = [] # --- NonIterative section ---# for group1 in propka_groups: for group2 in propka_groups: @@ -77,7 +83,7 @@ def add_determinants(group1, group2, distance, version): add_coulomb_determinants(group1, group2, distance, version) -def add_sidechain_determinants(group1, group2, version=None): +def add_sidechain_determinants(group1: Group, group2: Group, version: Version): """Add side-chain determinants and perturbations. NOTE - res_num1 > res_num2 @@ -236,6 +242,8 @@ def set_backbone_determinants(titratable_groups, backbone_groups, version): get_smallest_distance( backbone_interaction_atoms, titratable_group_interaction_atoms)) + assert backbone_atom is not None + assert titratable_atom is not None # get the parameters parameters = ( version.get_backbone_hydrogen_bond_parameters( diff --git a/propka/energy.py b/propka/energy.py index ae959ab..9f8c2a5 100644 --- a/propka/energy.py +++ b/propka/energy.py @@ -7,11 +7,15 @@ """ import math import logging -from typing import TYPE_CHECKING +from typing import TYPE_CHECKING, Optional, Sequence + +from propka.atom import Atom +from propka.parameters import Parameters if TYPE_CHECKING: from propka.conformation_container import ConformationContainer from propka.group import Group + from propka.version import Version from propka.calculations import squared_distance, get_smallest_distance @@ -89,7 +93,7 @@ def calculate_scale_factor(parameters, weight: float) -> float: return scale_factor -def calculate_weight(parameters, num_volume: int) -> float: +def calculate_weight(parameters: Parameters, num_volume: float) -> float: """Calculate the atom-based desolvation weight. TODO - figure out why a similar function exists in version.py @@ -109,7 +113,7 @@ def calculate_weight(parameters, num_volume: int) -> float: return weight -def calculate_pair_weight(parameters, num_volume1: int, num_volume2: int) -> float: +def calculate_pair_weight(parameters: Parameters, num_volume1: int, num_volume2: int) -> float: """Calculate the atom-pair based desolvation weight. Args: @@ -148,7 +152,11 @@ def hydrogen_bond_energy(dist, dpka_max: float, cutoffs, f_angle=1.0) -> float: return abs(dpka) -def angle_distance_factors(atom1=None, atom2=None, atom3=None, center=None): +def angle_distance_factors( + atom1: Optional[Atom] = None, + atom2: Atom = None, # type: ignore[assignment] + atom3: Atom = None, # type: ignore[assignment] + center: Optional[Sequence[float]] = None): """Calculate distance and angle factors for three atoms for backbone interactions. @@ -182,6 +190,7 @@ def angle_distance_factors(atom1=None, atom2=None, atom3=None, center=None): dy_32 = dy_32/dist_23 dz_32 = dz_32/dist_23 if atom1 is None: + assert center is not None dx_21 = center[0] - atom2.x dy_21 = center[1] - atom2.y dz_21 = center[2] - atom2.z @@ -197,7 +206,7 @@ def angle_distance_factors(atom1=None, atom2=None, atom3=None, center=None): return dist_12, f_angle, dist_23 -def hydrogen_bond_interaction(group1, group2, version): +def hydrogen_bond_interaction(group1: "Group", group2: "Group", version: "Version"): """Calculate energy for hydrogen bond interactions between two groups. Args: @@ -213,7 +222,7 @@ def hydrogen_bond_interaction(group1, group2, version): [closest_atom1, dist, closest_atom2] = get_smallest_distance( atoms1, atoms2 ) - if None in [closest_atom1, closest_atom2]: + if closest_atom1 is None or closest_atom2 is None: _LOGGER.warning( 'Side chain interaction failed for {0:s} and {1:s}'.format( group1.label, group2.label)) @@ -297,7 +306,7 @@ def electrostatic_interaction(group1, group2, dist, version): return value -def check_coulomb_pair(parameters, group1, group2, dist): +def check_coulomb_pair(parameters: Parameters, group1: "Group", group2: "Group", dist: float) -> bool: """Checks if this Coulomb interaction should be done. NOTE - this is a propka2.0 hack diff --git a/propka/group.py b/propka/group.py index 65e09f4..4ff702c 100644 --- a/propka/group.py +++ b/propka/group.py @@ -14,6 +14,7 @@ from typing import cast, Dict, Iterable, List, NoReturn, Optional import propka.ligand +from propka.parameters import Parameters import propka.protonate from propka.atom import Atom from propka.ligand_pka_values import LigandPkaValues @@ -90,7 +91,7 @@ def __init__(self, atom: Atom): self.y = 0.0 self.z = 0.0 self.charge = 0 - self.parameters = None + self.parameters: Optional[Parameters] = None self.exclude_cys_from_results = False self.interaction_atoms_for_acids: List[Atom] = [] self.interaction_atoms_for_bases: List[Atom] = [] @@ -167,6 +168,7 @@ def share_determinants(self, others: Iterable["Group"]) -> None: Args: others: list of other groups """ + raise NotImplementedError("unused") # for each determinant type for other in others: if other == self: @@ -319,6 +321,7 @@ def clone(self): def setup(self): """Set up a group.""" + assert self.parameters is not None # set the charges if self.type in self.parameters.charge.keys(): self.charge = self.parameters.charge[self.type] @@ -402,7 +405,7 @@ def set_interaction_atoms(self, interaction_atoms_for_acids: List[Atom], ' {0:s}'.format( str(self.interaction_atoms_for_bases[i]))) - def get_interaction_atoms(self, interacting_group) -> List[Atom]: + def get_interaction_atoms(self, interacting_group: "Group") -> List[Atom]: """Get atoms involved in interaction with other group. Args: @@ -591,7 +594,7 @@ def calculate_folding_energy(self, parameters, ph=None, reference=None): ddg = ddg_neutral + ddg_low return ddg - def calculate_charge(self, _, ph=7.0, state='folded'): + def calculate_charge(self, _, ph: float = 7.0, state: str = 'folded') -> float: """Calculate the charge of the specified state at the specified pH. Args: @@ -609,7 +612,7 @@ def calculate_charge(self, _, ph=7.0, state='folded'): charge = self.charge*(conc_ratio/(1.0+conc_ratio)) return charge - def use_in_calculations(self): + def use_in_calculations(self) -> bool: """Indicate whether group should be included in results report. If --titrate_only option is specified, only residues that are @@ -1218,7 +1221,7 @@ def __init__(self, atom): self.model_pka_set = True -def is_group(parameters, atom: Atom) -> Optional[Group]: +def is_group(parameters: Parameters, atom: Atom) -> Optional[Group]: """Identify whether the atom belongs to a group. Args: @@ -1227,7 +1230,7 @@ def is_group(parameters, atom: Atom) -> Optional[Group]: Returns: group for atom or None """ - atom.groups_extracted = 1 + atom.groups_extracted = True # check if this atom belongs to a protein group protein_group = is_protein_group(parameters, atom) if protein_group: @@ -1245,7 +1248,7 @@ def is_group(parameters, atom: Atom) -> Optional[Group]: ligand_group = is_ligand_group_by_groups(parameters, atom) else: raise Exception( - 'Unknown ligand typing method \'{0.s}\''.format( + 'Unknown ligand typing method \'{0:s}\''.format( parameters.ligand_typing)) if ligand_group: return ligand_group @@ -1368,7 +1371,7 @@ def is_ligand_group_by_groups(_, atom: Atom) -> Optional[Group]: return None -def is_ligand_group_by_marvin_pkas(parameters, atom: Atom) -> Optional[Group]: +def is_ligand_group_by_marvin_pkas(parameters: Parameters, atom: Atom) -> Optional[Group]: """Identify whether the atom belongs to a ligand group by calculating 'Marvin pKas'. @@ -1383,6 +1386,7 @@ def is_ligand_group_by_marvin_pkas(parameters, atom: Atom) -> Optional[Group]: # calculate Marvin ligand pkas for this conformation container # if not already done # TODO - double-check testing coverage of these functions. + assert atom.molecular_container is not None assert atom.conformation_container is not None if not atom.conformation_container.marvin_pkas_calculated: lpka = LigandPkaValues(parameters) @@ -1400,6 +1404,7 @@ def is_ligand_group_by_marvin_pkas(parameters, atom: Atom) -> Optional[Group]: if not o == atom][0] atom.marvin_pka = other_oxygen.marvin_pka return TitratableLigandGroup(atom) + raise NotImplementedError("hydrogen_bonds") if atom.element in parameters.hydrogen_bonds.elements: return NonTitratableLigandGroup(atom) return None diff --git a/propka/hydrogens.py b/propka/hydrogens.py index b7034ef..72fed3c 100644 --- a/propka/hydrogens.py +++ b/propka/hydrogens.py @@ -122,6 +122,12 @@ def add_arg_hydrogen(residue: List[Atom]) -> List[Atom]: Returns: list of hydrogen atoms """ + cd_atom: Optional[Atom] = None + cz_atom: Optional[Atom] = None + ne_atom: Optional[Atom] = None + nh1_atom: Optional[Atom] = None + nh2_atom: Optional[Atom] = None + for atom in residue: if atom.name == "CD": cd_atom = atom @@ -133,6 +139,11 @@ def add_arg_hydrogen(residue: List[Atom]) -> List[Atom]: nh1_atom = atom elif atom.name == "NH2": nh2_atom = atom + + if (cd_atom is None or cz_atom is None or ne_atom is None or nh1_atom is None + or nh2_atom is None): + raise ValueError("Unable to find all atoms") + h1_atom = protonate_sp2(cd_atom, ne_atom, cz_atom) h1_atom.name = "HE" h2_atom = protonate_direction(nh1_atom, ne_atom, cz_atom) @@ -152,6 +163,12 @@ def add_his_hydrogen(residue: List[Atom]) -> None: Args: residue: histidine residue to protonate """ + cg_atom: Optional[Atom] = None + nd_atom: Optional[Atom] = None + cd_atom: Optional[Atom] = None + ce_atom: Optional[Atom] = None + ne_atom: Optional[Atom] = None + for atom in residue: if atom.name == "CG": cg_atom = atom @@ -163,6 +180,11 @@ def add_his_hydrogen(residue: List[Atom]) -> None: ce_atom = atom elif atom.name == "NE2": ne_atom = atom + + if (cg_atom is None or nd_atom is None or cd_atom is None or ce_atom is None + or ne_atom is None): + raise ValueError("Unable to find all atoms") + hd_atom = protonate_sp2(cg_atom, nd_atom, ce_atom) hd_atom.name = "HND" he_atom = protonate_sp2(cd_atom, ne_atom, ce_atom) @@ -177,6 +199,7 @@ def add_trp_hydrogen(residue: List[Atom]) -> None: """ cd_atom = None ne_atom = None + ce_atom = None for atom in residue: if atom.name == "CD1": cd_atom = atom diff --git a/propka/input.py b/propka/input.py index e5658d8..6925d4c 100644 --- a/propka/input.py +++ b/propka/input.py @@ -9,8 +9,7 @@ Methods to read PROPKA input files (:func:`read_propka` and :func:`get_atom_lines_from_input`) have been removed. """ -import typing -from typing import Iterator, Tuple, Union +from typing import IO, ContextManager, Dict, Iterable, Iterator, Optional, Tuple import contextlib import io import zipfile @@ -19,19 +18,18 @@ from propka.atom import Atom from propka.conformation_container import ConformationContainer from propka.molecular_container import MolecularContainer +from propka.output import _PathArg, _PathLikeTypes, _TextIOSource from propka.parameters import Parameters -def open_file_for_reading( - input_file: typing.Union[str, Path, typing.TextIO] -) -> typing.ContextManager[typing.TextIO]: +def open_file_for_reading(input_file: _TextIOSource) -> ContextManager[IO[str]]: """Open file or file-like stream for reading. Args: input_file: path to file or file-like object. If file-like object, then will attempt seek(0). """ - if not isinstance(input_file, (str, Path)): + if not isinstance(input_file, _PathLikeTypes): input_file.seek(0) return contextlib.nullcontext(input_file) @@ -48,7 +46,11 @@ def open_file_for_reading( return contextlib.closing(open(input_file, 'rt')) -def read_molecule_file(filename: str, mol_container: MolecularContainer, stream=None) -> MolecularContainer: +def read_molecule_file( + filename: _PathArg, + mol_container: MolecularContainer, + stream: Optional[IO[str]] = None, +) -> MolecularContainer: """Read input file or stream (PDB or PROPKA) for a molecular container Args: @@ -96,11 +98,7 @@ def read_molecule_file(filename: str, mol_container: MolecularContainer, stream= input_path = Path(filename) mol_container.name = input_path.stem input_file_extension = input_path.suffix - - if stream is not None: - input_file = stream - else: - input_file = filename + input_file = filename if stream is None else stream if input_file_extension.lower() == '.pdb': # input is a pdb file. read in atoms and top up containers to make @@ -133,7 +131,7 @@ def read_molecule_file(filename: str, mol_container: MolecularContainer, stream= return mol_container -def read_parameter_file(input_file: Union[Path, str], parameters: Parameters) -> Parameters: +def read_parameter_file(input_file: _PathArg, parameters: Parameters) -> Parameters: """Read a parameter file. Args: @@ -161,8 +159,13 @@ def conformation_sorter(conf: str) -> int: return model*100+ord(altloc) -def get_atom_lines_from_pdb(pdb_file, ignore_residues=[], keep_protons=False, - tags=['ATOM ', 'HETATM'], chains=None) -> Iterator[Tuple[str, Atom]]: +def get_atom_lines_from_pdb( + pdb_file: _TextIOSource, + ignore_residues: Iterable[str] = (), + keep_protons: bool = False, + tags: Iterable[str] = ('ATOM ', 'HETATM'), + chains: Optional[Iterable[str]] = None, +) -> Iterator[Tuple[str, Atom]]: """Get atom lines from PDB file. Args: @@ -228,7 +231,8 @@ def get_atom_lines_from_pdb(pdb_file, ignore_residues=[], keep_protons=False, terminal = None -def read_pdb(pdb_file, parameters, molecule): +def read_pdb(pdb_file: _TextIOSource, parameters: Parameters, + molecule: MolecularContainer): """Parse a PDB file. Args: @@ -240,7 +244,7 @@ def read_pdb(pdb_file, parameters, molecule): 1. list of conformations 2. list of names """ - conformations = {} + conformations: Dict[str, ConformationContainer] = {} # read in all atoms in the file lines = get_atom_lines_from_pdb( pdb_file, ignore_residues=parameters.ignore_residues, @@ -253,4 +257,4 @@ def read_pdb(pdb_file, parameters, molecule): conformations[name].add_atom(atom) # make a sorted list of conformation names names = sorted(conformations.keys(), key=conformation_sorter) - return [conformations, names] + return conformations, names diff --git a/propka/iterative.py b/propka/iterative.py index 47b3595..2ada994 100644 --- a/propka/iterative.py +++ b/propka/iterative.py @@ -7,7 +7,10 @@ """ import logging +from typing import Dict, Iterable, List, Optional, Sequence, Tuple from propka.determinant import Determinant +from propka.group import Group +from propka.version import Version _LOGGER = logging.getLogger(__name__) @@ -16,9 +19,11 @@ # TODO - these are undocumented constants UNK_MIN_VALUE = 0.005 +Interaction = list -def add_to_determinant_list(group1, group2, distance, iterative_interactions, - version): + +def add_to_determinant_list(group1: Group, group2: Group, distance: float, + iterative_interactions: List[Interaction], version: Version): """Add iterative determinantes to the list. [[R1, R2], [side-chain, coulomb], [A1, A2]], ... @@ -48,7 +53,8 @@ def add_to_determinant_list(group1, group2, distance, iterative_interactions, iterative_interactions.append(interaction) -def add_iterative_acid_pair(object1, object2, interaction): +def add_iterative_acid_pair(object1: "Iterative", object2: "Iterative", + interaction: Interaction): """Add the Coulomb 'iterative' interaction (an acid pair). The higher pKa is raised with QQ+HB @@ -90,7 +96,8 @@ def add_iterative_acid_pair(object1, object2, interaction): annihilation[1] = -diff -def add_iterative_base_pair(object1, object2, interaction): +def add_iterative_base_pair(object1: "Iterative", object2: "Iterative", + interaction: Interaction): """Add the Coulomb 'iterative' interaction (a base pair). The lower pKa is lowered @@ -132,7 +139,8 @@ def add_iterative_base_pair(object1, object2, interaction): annihilation[1] = -diff -def add_iterative_ion_pair(object1, object2, interaction, version): +def add_iterative_ion_pair(object1: "Iterative", object2: "Iterative", + interaction: Interaction, version: Version): """Add the Coulomb 'iterative' interaction (an acid-base pair) the pKa of the acid is lowered & the pKa of the base is raised @@ -194,7 +202,7 @@ def add_iterative_ion_pair(object1, object2, interaction, version): object2.determinants['sidechain'].append(interaction) -def add_determinants(iterative_interactions, version, _=None): +def add_determinants(iterative_interactions: List[Interaction], version: Version, _=None): """Add determinants iteratively. The iterative pKa scheme. Later it is all added in 'calculateTotalPKA' @@ -205,7 +213,7 @@ def add_determinants(iterative_interactions, version, _=None): _: options object """ # --- setup --- - iteratives = [] + iteratives: List[Iterative] = [] done_group = [] # create iterative objects with references to their real group counterparts for interaction in iterative_interactions: @@ -270,6 +278,7 @@ def add_determinants(iterative_interactions, version, _=None): # reset pka_old & storing pka_new in pka_iter for itres in iteratives: + assert itres.pka_new is not None itres.pka_old = itres.pka_new itres.pka_iter.append(itres.pka_new) @@ -295,14 +304,17 @@ def add_determinants(iterative_interactions, version, _=None): for itres in iteratives: for type_ in ['sidechain', 'backbone', 'coulomb']: for interaction in itres.determinants[type_]: - value = interaction[1] + value: float = interaction[1] if value > UNK_MIN_VALUE or value < -UNK_MIN_VALUE: group = interaction[0] new_det = Determinant(group, value) itres.group.determinants[type_].append(new_det) -def find_iterative(pair, iteratives): +def find_iterative( + pair: Sequence[Group], + iteratives: Iterable["Iterative"], +) -> Tuple["Iterative", "Iterative"]: """Find the 'iteratives' that correspond to the groups in 'pair'. Args: @@ -312,11 +324,15 @@ def find_iterative(pair, iteratives): 1. first matched iterative 2. second matched iterative """ + iterative0: Optional[Iterative] = None + iterative1: Optional[Iterative] = None for iterative in iteratives: if iterative.group == pair[0]: iterative0 = iterative elif iterative.group == pair[1]: iterative1 = iterative + if iterative0 is None or iterative1 is None: + raise LookupError("iteratives not found") return iterative0, iterative1 @@ -327,7 +343,7 @@ class Iterative: after the iterations are finished. """ - def __init__(self, group): + def __init__(self, group: Group): """Initialize object with group. Args: @@ -337,11 +353,15 @@ def __init__(self, group): self.atom = group.atom self.res_name = group.residue_type self.q = group.charge - self.pka_old = None - self.pka_new = None - self.pka_iter = [] + self.pka_old: Optional[float] = None + self.pka_new: Optional[float] = None + self.pka_iter: List[float] = [] self.pka_noniterative = 0.00 - self.determinants = {'sidechain': [], 'backbone': [], 'coulomb': []} + self.determinants: Dict[str, list] = { + 'sidechain': [], + 'backbone': [], + 'coulomb': [] + } self.group = group self.converged = True # Calculate the Non-Iterative part of pKa from the group object @@ -368,8 +388,9 @@ def __init__(self, group): self.pka_noniterative += coulomb self.pka_old = self.pka_noniterative - def __eq__(self, other): + def __eq__(self, other) -> bool: """Needed to use objects in sets.""" + assert isinstance(other, (Iterative, Group)), type(other) if self.atom.type == 'atom': # In case of protein atoms we trust the labels return self.label == other.label diff --git a/propka/lib.py b/propka/lib.py index 4519094..17dfc86 100644 --- a/propka/lib.py +++ b/propka/lib.py @@ -5,11 +5,19 @@ Implements many of the main functions used to call PROPKA. """ -import sys import logging import argparse from pathlib import Path +from typing import Iterable, Iterator, List, TYPE_CHECKING, NoReturn, Optional, Tuple, TypeVar +if TYPE_CHECKING: + from propka.atom import Atom + + +T = TypeVar("T") +Number = TypeVar("Number", int, float) + +_T_RESIDUE_TUPLE = Tuple[str, int, str] _LOGGER = logging.getLogger(__name__) @@ -20,6 +28,29 @@ 'SER': 6, 'THR': 7, 'TRP': 14, 'TYR': 12, 'VAL': 7} +class Options: + # Note: All the "NoReturn" members appear to be unused + alignment: NoReturn # Optional[List[str]] + chains: Optional[List[str]] + display_coupled_residues: bool = False + filenames: List[str] # List[Path]? + grid: Tuple[float, float, float] = (0.0, 14.0, 0.1) + input_pdb: str # Path? + keep_protons: bool = False + log_level: str = 'INFO' + mutations: NoReturn # Optional[List[str]] + mutator: NoReturn # Optional[str] # alignment/scwrl/jackal + mutator_options: NoReturn # Optional[List[str]] + pH: NoReturn # float = 7.0 + parameters: Path + protonate_all: bool = False + reference: NoReturn # str = 'neutral' + reuse_ligand_mol2_file: bool = False # only used by unused function + thermophiles: NoReturn # Optional[List[str]] + titrate_only: Optional[List[_T_RESIDUE_TUPLE]] + window: Tuple[float, float, float] = (0.0, 14.0, 1.0) + + def protein_precheck(conformations, names): """Check protein for correct number of atoms, etc. @@ -73,7 +104,7 @@ def resid_from_atom(atom): atom.res_num, atom.chain_id, atom.icode) -def split_atoms_into_molecules(atoms): +def split_atoms_into_molecules(atoms: List["Atom"]): """Maps atoms into molecules. Args: @@ -81,14 +112,14 @@ def split_atoms_into_molecules(atoms): Returns: list of molecules """ - molecules = [] + molecules: List[List["Atom"]] = [] while len(atoms) > 0: initial_atom = atoms.pop() molecules.append(make_molecule(initial_atom, atoms)) return molecules -def make_molecule(atom, atoms): +def make_molecule(atom: "Atom", atoms: List["Atom"]): """Make a molecule from atoms. Args: @@ -106,11 +137,11 @@ def make_molecule(atom, atoms): return res_atoms -def make_grid(min_, max_, step): +def make_grid(min_: Number, max_: Number, step: Number) -> Iterator[Number]: """Make a grid across the specified tange. - TODO - figure out if this duplicates existing generators like `range` or - numpy function. + Like range() for integers or numpy.arange() for floats, except that `max_` + is not excluded from the range. Args: min_: minimum value of grid @@ -123,7 +154,7 @@ def make_grid(min_, max_, step): x += step -def generate_combinations(interactions): +def generate_combinations(interactions: Iterable[T]) -> List[List[T]]: """Generate combinations of interactions. Args: @@ -131,14 +162,14 @@ def generate_combinations(interactions): Returns: list of combinations """ - res = [[]] + res: List[List[T]] = [[]] for interaction in interactions: res = make_combination(res, interaction) res.remove([]) return res -def make_combination(combis, interaction): +def make_combination(combis: List[List[T]], interaction: T) -> List[List[T]]: """Make a specific set of combinations. Args: @@ -154,7 +185,7 @@ def make_combination(combis, interaction): return res -def parse_res_string(res_str): +def parse_res_string(res_str: str) -> _T_RESIDUE_TUPLE: """Parse a residue string. Args: @@ -183,6 +214,16 @@ def parse_res_string(res_str): return chain, resnum, inscode +def parse_res_list(titrate_only: str): + res_list: List[_T_RESIDUE_TUPLE] = [] + for res_str in titrate_only.split(','): + try: + res_list.append(parse_res_string(res_str)) + except ValueError as ex: + raise argparse.ArgumentTypeError(f'{ex}: "{res_str:s}"') + return res_list + + def build_parser(parser=None): """Build an argument parser for PROPKA. @@ -227,6 +268,7 @@ def build_parser(parser=None): '" " for chains without ID [all]')) group.add_argument( "-i", "--titrate_only", dest="titrate_only", + type=parse_res_list, help=('Treat only the specified residues as titratable. Value should ' 'be a comma-separated list of "chain:resnum" values; for ' 'example: -i "A:10,A:11"')) @@ -246,8 +288,8 @@ def build_parser(parser=None): "--version", action="version", version=f"%(prog)s {propka.__version__}") group.add_argument( "-p", "--parameters", dest="parameters", - default=str(Path(__file__).parent / "propka.cfg"), - help="set the parameter file [{default:s}]") + type=Path, default=Path(__file__).parent / "propka.cfg", + help="set the parameter file") try: group.add_argument( "--log-level", @@ -302,7 +344,7 @@ def build_parser(parser=None): return parser -def loadOptions(args=None): +def loadOptions(args=None) -> Options: """ Load the arguments parser with options. Note that verbosity is set as soon as this function is invoked. @@ -315,22 +357,10 @@ def loadOptions(args=None): # loading the parser parser = build_parser() # parsing and returning options and arguments - options = parser.parse_args(args) + options = parser.parse_args(args, namespace=Options()) # adding specified filenames to arguments options.filenames.append(options.input_pdb) - # Convert titrate_only string to a list of (chain, resnum) items: - if options.titrate_only is not None: - res_list = [] - for res_str in options.titrate_only.split(','): - try: - chain, resnum, inscode = parse_res_string(res_str) - except ValueError: - _LOGGER.critical( - 'Invalid residue string: "{0:s}"'.format(res_str)) - sys.exit(1) - res_list.append((chain, resnum, inscode)) - options.titrate_only = res_list # Set the no-print variable level = getattr(logging, options.log_level) _LOGGER.setLevel(level) diff --git a/propka/ligand.py b/propka/ligand.py index 618f105..88f389b 100644 --- a/propka/ligand.py +++ b/propka/ligand.py @@ -319,11 +319,11 @@ def are_atoms_planar(atoms): return False vec1 = Vector(atom1=atoms[0], atom2=atoms[1]) vec2 = Vector(atom1=atoms[0], atom2=atoms[2]) - norm = (vec1**vec2).rescale(1.0) + norm = vec1.cross(vec2).rescale(1.0) margin = PLANARITY_MARGIN for atom in atoms[3:]: vec = Vector(atom1=atoms[0], atom2=atom).rescale(1.0) - if abs(vec*norm) > margin: + if abs(vec.dot(norm)) > margin: return False return True diff --git a/propka/ligand_pka_values.py b/propka/ligand_pka_values.py index f161d88..c36f4d2 100644 --- a/propka/ligand_pka_values.py +++ b/propka/ligand_pka_values.py @@ -11,11 +11,18 @@ """ import logging import os +import shutil import subprocess import sys +import warnings +from typing import TYPE_CHECKING, NoReturn + from propka.output import write_mol2_for_atoms from propka.lib import split_atoms_into_molecules +from propka.parameters import Parameters +if TYPE_CHECKING: + from propka.molecular_container import MolecularContainer _LOGGER = logging.getLogger(__name__) @@ -23,7 +30,7 @@ class LigandPkaValues: """Ligand pKa value class.""" - def __init__(self, parameters): + def __init__(self, parameters: Parameters): """Initialize object with parameters. Args: @@ -46,20 +53,17 @@ def find_in_path(program): Returns: location of program """ - path = os.environ.get('PATH').split(os.pathsep) - locs = [ - i for i in filter(lambda loc: os.access(loc, os.F_OK), - map(lambda dir: os.path.join(dir, program), - path))] - if len(locs) == 0: + loc = shutil.which(program) + if loc is None: str_ = "'Error: Could not find {0:s}.".format(program) str_ += ' Please make sure that it is found in the path.' _LOGGER.info(str_) sys.exit(-1) - return locs[0] + return loc def get_marvin_pkas_for_pdb_file( - self, molecule, parameters, num_pkas=10, min_ph=-10, max_ph=20): + self, molecule: "MolecularContainer", parameters: NoReturn, + num_pkas=10, min_ph=-10.0, max_ph=20.0): """Use Marvin executables to get pKas for a PDB file. Args: @@ -69,11 +73,12 @@ def get_marvin_pkas_for_pdb_file( min_ph: minimum pH value max_ph: maximum pH value """ + warnings.warn("unused and untested by propka") self.get_marvin_pkas_for_molecular_container( molecule, num_pkas=num_pkas, min_ph=min_ph, max_ph=max_ph) - def get_marvin_pkas_for_molecular_container(self, molecule, num_pkas=10, - min_ph=-10, max_ph=20): + def get_marvin_pkas_for_molecular_container(self, molecule: "MolecularContainer", num_pkas=10, + min_ph=-10.0, max_ph=20.0): """Use Marvin executables to calculate pKas for a molecular container. Args: @@ -91,8 +96,8 @@ def get_marvin_pkas_for_molecular_container(self, molecule, num_pkas=10, def get_marvin_pkas_for_conformation_container(self, conformation, name='temp', reuse=False, - num_pkas=10, min_ph=-10, - max_ph=20): + num_pkas=10, min_ph=-10.0, + max_ph=20.0): """Use Marvin executables to calculate pKas for a conformation container. Args: @@ -109,7 +114,7 @@ def get_marvin_pkas_for_conformation_container(self, conformation, num_pkas=num_pkas, min_ph=min_ph, max_ph=max_ph) def get_marvin_pkas_for_atoms(self, atoms, name='temp', reuse=False, - num_pkas=10, min_ph=-10, max_ph=20): + num_pkas=10, min_ph=-10.0, max_ph=20.0): """Use Marvin executables to calculate pKas for a list of atoms. Args: @@ -129,8 +134,8 @@ def get_marvin_pkas_for_atoms(self, atoms, name='temp', reuse=False, min_ph=min_ph, max_ph=max_ph) def get_marvin_pkas_for_molecule(self, atoms, filename='__tmp_ligand.mol2', - reuse=False, num_pkas=10, min_ph=-10, - max_ph=20): + reuse=False, num_pkas=10, min_ph=-10.0, + max_ph=20.0): """Use Marvin executables to calculate pKas for a molecule. Args: diff --git a/propka/molecular_container.py b/propka/molecular_container.py index 2f7f892..394ae6f 100644 --- a/propka/molecular_container.py +++ b/propka/molecular_container.py @@ -7,11 +7,12 @@ import logging import os from typing import Dict, List, Optional, Tuple +from propka.parameters import Parameters import propka.version from propka.output import write_pka, print_header, print_result from propka.conformation_container import ConformationContainer -from propka.lib import make_grid +from propka.lib import make_grid, Options _LOGGER = logging.getLogger(__name__) @@ -35,7 +36,7 @@ class MolecularContainer: name: Optional[str] version: propka.version.Version - def __init__(self, parameters, options=None) -> None: + def __init__(self, parameters: Parameters, options: Options) -> None: """Initialize molecular container. Args: @@ -156,7 +157,7 @@ def write_pka(self, filename=None, reference="neutral", conformation='AVR', reference=reference) def get_folding_profile(self, conformation='AVR', reference="neutral", - grid=[0., 14., 0.1]): + grid: Tuple[float, float, float] = (0., 14., 0.1)): """Get a folding profile. Args: @@ -173,25 +174,25 @@ def get_folding_profile(self, conformation='AVR', reference="neutral", 4. stability_range """ # calculate stability profile - profile = [] + profile: List[Tuple[float, float]] = [] for ph in make_grid(*grid): conf = self.conformations[conformation] ddg = conf.calculate_folding_energy(ph=ph, reference=reference) - profile.append([ph, ddg]) + profile.append((ph, ddg)) # find optimum - opt = [None, 1e6] + opt: Tuple[Optional[float], float] = (None, 1e6) for point in profile: opt = min(opt, point, key=lambda v: v[1]) # find values within 80 % of optimum - range_80pct = [None, None] + range_80pct: Tuple[Optional[float], Optional[float]] = (None, None) values_within_80pct = [p[0] for p in profile if p[1] < 0.8*opt[1]] if len(values_within_80pct) > 0: - range_80pct = [min(values_within_80pct), max(values_within_80pct)] + range_80pct = (min(values_within_80pct), max(values_within_80pct)) # find stability range - stability_range = [None, None] + stability_range: Tuple[Optional[float], Optional[float]] = (None, None) stable_values = [p[0] for p in profile if p[1] < 0.0] if len(stable_values) > 0: - stability_range = [min(stable_values), max(stable_values)] + stability_range = (min(stable_values), max(stable_values)) return profile, opt, range_80pct, stability_range def get_charge_profile(self, conformation: str = 'AVR', grid=[0., 14., .1]): diff --git a/propka/output.py b/propka/output.py index db5c311..d9d8617 100644 --- a/propka/output.py +++ b/propka/output.py @@ -12,13 +12,29 @@ import logging from datetime import date from decimal import Decimal +from os import PathLike +from pathlib import Path +from typing import IO, AnyStr, List, NoReturn, Optional, Union, TYPE_CHECKING +import warnings + +from .parameters import Parameters from . import __version__ +if TYPE_CHECKING: + from .atom import Atom + from .conformation_container import ConformationContainer + from .molecular_container import MolecularContainer + +# https://docs.python.org/3/glossary.html#term-path-like-object +_PathLikeTypes = (PathLike, str) +_PathArg = Union[PathLike, str] +_IOSource = Union[IO[AnyStr], PathLike, str] +_TextIOSource = _IOSource[str] _LOGGER = logging.getLogger(__name__) -def open_file_for_writing(input_file): +def open_file_for_writing(input_file: _TextIOSource) -> IO[str]: """Open file or file-like stream for writing. TODO - convert this to a context manager. @@ -27,17 +43,13 @@ def open_file_for_writing(input_file): input_file: path to file or file-like object. If file-like object, then will attempt to get file mode. """ - try: - if not input_file.writable(): - raise IOError("File/stream not open for writing") - return input_file - except AttributeError: - pass - try: - file_ = open(input_file, 'wt') - except FileNotFoundError: - raise Exception('Could not open {0:s}'.format(input_file)) - return file_ + if isinstance(input_file, _PathLikeTypes): + return open(input_file, 'wt') + + if not input_file.writable(): + raise IOError("File/stream not open for writing") + + return input_file def write_file(filename, lines): @@ -47,6 +59,7 @@ def write_file(filename, lines): filename: name of file lines: lines to write to file """ + warnings.warn("unused and untested by propka") file_ = open_file_for_writing(filename) for line in lines: file_.write("{0:s}\n".format(line)) @@ -72,6 +85,7 @@ def write_pdb_for_protein( include_hydrogens: Boolean indicating whether to include hydrogens options: options object """ + raise NotImplementedError("unused") if pdbfile is None: # opening file if not given if filename is None: @@ -100,17 +114,22 @@ def write_pdb_for_protein( pdbfile.close() -def write_pdb_for_conformation(conformation, filename): +def write_pdb_for_conformation(conformation: "ConformationContainer", + filename: _PathArg): """Write PDB conformation to a file. Args: conformation: conformation container filename: filename for output """ + warnings.warn("unused and untested by propka") write_pdb_for_atoms(conformation.atoms, filename) -def write_pka(protein, parameters, filename=None, conformation='1A', +def write_pka(protein: "MolecularContainer", + parameters: Parameters, + filename: Optional[_PathArg] = None, + conformation='1A', reference="neutral", _="folding", verbose=False, __=None): """Write the pKa-file based on the given protein. @@ -128,8 +147,6 @@ def write_pka(protein, parameters, filename=None, conformation='1A', verbose = True if filename is None: filename = "{0:s}.pka".format(protein.name) - # TODO - this would be much better with a context manager - file_ = open(filename, 'w') if verbose: _LOGGER.info("Writing {0:s}".format(filename)) # writing propka header @@ -148,11 +165,10 @@ def write_pka(protein, parameters, filename=None, conformation='1A', # printing Protein Charge Profile str_ += get_charge_profile_section(protein, conformation=conformation) # now, writing the pka text to file - file_.write(str_) - file_.close() + Path(filename).write_text(str_, encoding="utf-8") -def print_tm_profile(protein, reference="neutral", window=[0., 14., 1.], +def print_tm_profile(protein: NoReturn, reference="neutral", window=[0., 14., 1.], __=[0., 0.], tms=None, ref=None, _=False, options=None): """Print Tm profile. @@ -169,6 +185,7 @@ def print_tm_profile(protein, reference="neutral", window=[0., 14., 1.], _: Boolean for verbosity options: options object """ + raise NotImplementedError("unused") profile = protein.getTmProfile( reference=reference, grid=[0., 14., 0.1], tms=tms, ref=ref, options=options) @@ -184,7 +201,7 @@ def print_tm_profile(protein, reference="neutral", window=[0., 14., 1.], _LOGGER.info(str_) -def print_result(protein, conformation, parameters): +def print_result(protein: "MolecularContainer", conformation: str, parameters: Parameters): """Prints all resulting output from determinants and down. Args: @@ -195,7 +212,7 @@ def print_result(protein, conformation, parameters): print_pka_section(protein, conformation, parameters) -def print_pka_section(protein, conformation, parameters): +def print_pka_section(protein: "MolecularContainer", conformation: str, parameters: Parameters): """Prints out pKa section of results. Args: @@ -210,7 +227,7 @@ def print_pka_section(protein, conformation, parameters): _LOGGER.info("pKa summary:\n%s", str_) -def get_determinant_section(protein, conformation, parameters): +def get_determinant_section(protein: "MolecularContainer", conformation: str, parameters: Parameters): """Returns string with determinant section of results. Args: @@ -242,7 +259,8 @@ def get_determinant_section(protein, conformation, parameters): return str_ -def get_summary_section(protein, conformation, parameters): +def get_summary_section(protein: "MolecularContainer", conformation: str, + parameters: Parameters): """Returns string with summary section of the results. Args: @@ -264,7 +282,8 @@ def get_summary_section(protein, conformation, parameters): def get_folding_profile_section( - protein, conformation='AVR', direction="folding", reference="neutral", + protein: "MolecularContainer", + conformation='AVR', direction="folding", reference="neutral", window=[0., 14., 1.0], _=False, __=None): """Returns string with the folding profile section of the results. @@ -358,11 +377,12 @@ def write_jackal_scap_file(mutation_data=None, filename="1xxx_scap.list", TODO - figure out what this is """ + raise NotImplementedError("unused") with open(filename, 'w') as file_: for chain_id, _, res_num, code2 in mutation_data: str_ = "{chain:s}, {num:d}, {code:s}\n".format( chain=chain_id, num=res_num, code=code2) - file_.write(str_) + file_.write(str_) def write_scwrl_sequence_file(sequence, filename="x-ray.seq", _=None): @@ -370,6 +390,7 @@ def write_scwrl_sequence_file(sequence, filename="x-ray.seq", _=None): TODO - figure out what this is """ + warnings.warn("unused and untested by propka") with open(filename, 'w') as file_: start = 0 while len(sequence[start:]) > 60: @@ -533,7 +554,7 @@ def make_interaction_map(name, list_, interaction): return res -def write_pdb_for_atoms(atoms, filename, make_conect_section=False): +def write_pdb_for_atoms(atoms: List["Atom"], filename: _PathArg, make_conect_section=False): """Write out PDB file for atoms. Args: diff --git a/propka/parameters.py b/propka/parameters.py index aa9abcc..654b3fa 100644 --- a/propka/parameters.py +++ b/propka/parameters.py @@ -11,67 +11,128 @@ """ import logging +from dataclasses import dataclass, field +from typing import Dict, List + +try: + # New in version 3.10, deprecated since version 3.12 + from typing import TypeAlias +except ImportError: + TypeAlias = "TypeAlias" # type: ignore _LOGGER = logging.getLogger(__name__) -#: matrices -MATRICES = ['interaction_matrix'] -#: pari-wise matrices -PAIR_WISE_MATRICES = ['sidechain_cutoffs'] -#: :class:`dict` containing numbers -NUMBER_DICTIONARIES = [ - 'VanDerWaalsVolume', 'charge', 'model_pkas', 'ions', 'valence_electrons', - 'custom_model_pkas'] -#: :class:`dict` containing lists -LIST_DICTIONARIES = ['backbone_NH_hydrogen_bond', 'backbone_CO_hydrogen_bond'] -#: :class:`dict` containing strings -STRING_DICTIONARIES = ['protein_group_mapping'] -#: :class:`list` containing strings -STRING_LISTS = [ - 'ignore_residues', 'angular_dependent_sidechain_interactions', - 'acid_list', 'base_list', 'exclude_sidechain_interactions', - 'backbone_reorganisation_list', 'write_out_order'] -#: distances (:class:`float`) -DISTANCES = ['desolv_cutoff', 'buried_cutoff', 'coulomb_cutoff1', - 'coulomb_cutoff2'] -#: other parameters -PARAMETERS = [ - 'Nmin', 'Nmax', 'desolvationSurfaceScalingFactor', 'desolvationPrefactor', - 'desolvationAllowance', 'coulomb_diel', 'COO_HIS_exception', - 'OCO_HIS_exception', 'CYS_HIS_exception', 'CYS_CYS_exception', - 'min_ligand_model_pka', 'max_ligand_model_pka', - 'include_H_in_interactions', 'coupling_max_number_of_bonds', - 'min_bond_distance_for_hydrogen_bonds', 'coupling_penalty', - 'shared_determinants', 'common_charge_centre', 'hide_penalised_group', - 'remove_penalised_group', 'max_intrinsic_pka_diff', - 'min_interaction_energy', 'max_free_energy_diff', 'min_swap_pka_shift', - 'min_pka', 'max_pka', 'sidechain_interaction'] -# :class:`str` parameters -STRINGS = ['version', 'output_file_tag', 'ligand_typing', 'pH', 'reference'] +class squared_property: + + def __set_name__(self, owner, name: str): + assert name.endswith("_squared") + self._name_not_squared = name[:-len("_squared")] # removesuffix() + + def __get__(self, instance, owner=None) -> float: + if instance is None: + return self # type: ignore[return-value] + return getattr(instance, self._name_not_squared)**2 + + def __set__(self, instance, value: float): + setattr(instance, self._name_not_squared, value**0.5) + + +_T_MATRIX: TypeAlias = "InteractionMatrix" +_T_PAIR_WISE_MATRIX: TypeAlias = "PairwiseMatrix" +_T_NUMBER_DICTIONARY = Dict[str, float] +_T_LIST_DICTIONARY = Dict[str, list] +_T_STRING_DICTIONARY = Dict[str, str] +_T_STRING_LIST = List[str] +_T_STRING = str +_T_BOOL = bool +@dataclass class Parameters: """PROPKA parameter class.""" - def __init__(self): - """Initialize parameter class. - - Args: - parameter_file: file with parameters - """ - # TODO - need to define all members explicitly - self.model_pkas = {} - self.interaction_matrix = InteractionMatrix("interaction_matrix") - self.sidechain_cutoffs = None - # TODO - it would be nice to rename these; they're defined everywhere - self.COO_HIS_exception = None - self.OCO_HIS_exception = None - self.CYS_HIS_exception = None - self.CYS_CYS_exception = None - # These functions set up remaining data structures implicitly - self.set_up_data_structures() + # MATRICES + interaction_matrix: _T_MATRIX = field( + default_factory=lambda: InteractionMatrix("interaction_matrix")) + + # PAIR_WISE_MATRICES + sidechain_cutoffs: _T_PAIR_WISE_MATRIX = field( + default_factory=lambda: PairwiseMatrix("sidechain_cutoffs")) + + # NUMBER_DICTIONARIES + VanDerWaalsVolume: _T_NUMBER_DICTIONARY = field(default_factory=dict) + charge: _T_NUMBER_DICTIONARY = field(default_factory=dict) + model_pkas: _T_NUMBER_DICTIONARY = field(default_factory=dict) + ions: _T_NUMBER_DICTIONARY = field(default_factory=dict) + valence_electrons: _T_NUMBER_DICTIONARY = field(default_factory=dict) + custom_model_pkas: _T_NUMBER_DICTIONARY = field(default_factory=dict) + + # LIST_DICTIONARIES + backbone_NH_hydrogen_bond: _T_LIST_DICTIONARY = field(default_factory=dict) + backbone_CO_hydrogen_bond: _T_LIST_DICTIONARY = field(default_factory=dict) + + # STRING_DICTIONARIES + protein_group_mapping: _T_STRING_DICTIONARY = field(default_factory=dict) + + # STRING_LISTS + ignore_residues: _T_STRING_LIST = field(default_factory=list) + angular_dependent_sidechain_interactions: _T_STRING_LIST = field(default_factory=list) + acid_list: _T_STRING_LIST = field(default_factory=list) + base_list: _T_STRING_LIST = field(default_factory=list) + exclude_sidechain_interactions: _T_STRING_LIST = field(default_factory=list) + backbone_reorganisation_list: _T_STRING_LIST = field(default_factory=list) + write_out_order: _T_STRING_LIST = field(default_factory=list) + + # DISTANCES + desolv_cutoff: float = 20.0 + buried_cutoff: float = 15.0 + coulomb_cutoff1: float = 4.0 + coulomb_cutoff2: float = 10.0 + + # DISTANCES SQUARED + desolv_cutoff_squared = squared_property() + buried_cutoff_squared = squared_property() + coulomb_cutoff1_squared = squared_property() + coulomb_cutoff2_squared = squared_property() + + # STRINGS + version: _T_STRING = "VersionA" + output_file_tag: _T_STRING = "" + ligand_typing: _T_STRING = "groups" + pH: _T_STRING = "variable" + reference: _T_STRING = "neutral" + + # PARAMETERS + Nmin: int = 280 + Nmax: int = 560 + desolvationSurfaceScalingFactor: float = 0.25 + desolvationPrefactor: float = -13.0 + desolvationAllowance: float = 0.0 + coulomb_diel: float = 80.0 + # TODO - it would be nice to rename these; they're defined everywhere + COO_HIS_exception: float = 1.60 + OCO_HIS_exception: float = 1.60 + CYS_HIS_exception: float = 1.60 + CYS_CYS_exception: float = 3.60 + min_ligand_model_pka: float = -10.0 + max_ligand_model_pka: float = 20.0 + # include_H_in_interactions: NoReturn = None + coupling_max_number_of_bonds: int = 3 + min_bond_distance_for_hydrogen_bonds: int = 4 + # coupling_penalty: NoReturn = None + shared_determinants: _T_BOOL = False + common_charge_centre: _T_BOOL = False + # hide_penalised_group: NoReturn = None + remove_penalised_group: _T_BOOL = True + max_intrinsic_pka_diff: float = 2.0 + min_interaction_energy: float = 0.5 + max_free_energy_diff: float = 1.0 + min_swap_pka_shift: float = 1.0 + min_pka: float = 0.0 + max_pka: float = 10.0 + sidechain_interaction: float = 0.85 def parse_line(self, line): """Parse parameter file line.""" @@ -84,22 +145,21 @@ def parse_line(self, line): if len(words) == 0: return # parse the words - if len(words) == 3 and words[0] in NUMBER_DICTIONARIES: + typeannotation = self.__annotations__.get(words[0]) + if typeannotation is _T_NUMBER_DICTIONARY: self.parse_to_number_dictionary(words) - elif len(words) == 2 and words[0] in STRING_LISTS: + elif typeannotation is _T_STRING_LIST: self.parse_to_string_list(words) - elif len(words) == 2 and words[0] in DISTANCES: - self.parse_distance(words) - elif len(words) == 2 and words[0] in PARAMETERS: - self.parse_parameter(words) - elif len(words) == 2 and words[0] in STRINGS: + elif typeannotation is _T_STRING: self.parse_string(words) - elif len(words) > 2 and words[0] in LIST_DICTIONARIES: + elif typeannotation is _T_LIST_DICTIONARY: self.parse_to_list_dictionary(words) - elif words[0] in MATRICES+PAIR_WISE_MATRICES: + elif typeannotation is _T_MATRIX or typeannotation is _T_PAIR_WISE_MATRIX: self.parse_to_matrix(words) - elif len(words) == 3 and words[0] in STRING_DICTIONARIES: + elif typeannotation is _T_STRING_DICTIONARY: self.parse_to_string_dictionary(words) + else: + self.parse_parameter(words) def parse_to_number_dictionary(self, words): """Parse field to number dictionary. @@ -107,6 +167,7 @@ def parse_to_number_dictionary(self, words): Args: words: strings to parse. """ + assert len(words) == 3, words dict_ = getattr(self, words[0]) key = words[1] value = words[2] @@ -118,17 +179,19 @@ def parse_to_string_dictionary(self, words): Args: words: strings to parse """ + assert len(words) == 3, words dict_ = getattr(self, words[0]) key = words[1] value = words[2] dict_[key] = value - def parse_to_list_dictionary(self, words): + def parse_to_list_dictionary(self, words: List[str]): """Parse field to list dictionary. Args: words: strings to parse. """ + assert len(words) > 2, words dict_ = getattr(self, words[0]) key = words[1] if key not in dict_: @@ -144,6 +207,7 @@ def parse_to_string_list(self, words): Args: words: strings to parse """ + assert len(words) == 2, words list_ = getattr(self, words[0]) value = words[1] list_.append(value) @@ -158,24 +222,13 @@ def parse_to_matrix(self, words): value = tuple(words[1:]) matrix.add(value) - def parse_distance(self, words): - """Parse field to distance. - - Args: - words: strings to parse - """ - value = float(words[1]) - setattr(self, words[0], value) - value_sq = value*value - attr = "{0:s}_squared".format(words[0]) - setattr(self, attr, value_sq) - def parse_parameter(self, words): """Parse field to parameters. Args: words: strings to parse """ + assert len(words) == 2, words value = float(words[1]) setattr(self, words[0], value) @@ -185,28 +238,9 @@ def parse_string(self, words): Args: words: strings to parse """ + assert len(words) == 2, words setattr(self, words[0], words[1]) - def set_up_data_structures(self): - """Set up internal data structures. - - TODO - it would be better to make these assignments explicit in - __init__. - """ - for key_word in (NUMBER_DICTIONARIES + LIST_DICTIONARIES - + STRING_DICTIONARIES): - setattr(self, key_word, {}) - for key_word in STRING_LISTS: - setattr(self, key_word, []) - for key_word in STRINGS: - setattr(self, key_word, "") - for key_word in MATRICES: - matrix = InteractionMatrix(key_word) - setattr(self, key_word, matrix) - for key_word in PAIR_WISE_MATRICES: - matrix = PairwiseMatrix(key_word) - setattr(self, key_word, matrix) - def print_interaction_parameters(self): """Print interaction parameters.""" _LOGGER.info('--------------- Model pKa values ----------------------') diff --git a/propka/propka.cfg b/propka/propka.cfg index c78ea43..93b09d4 100644 --- a/propka/propka.cfg +++ b/propka/propka.cfg @@ -107,7 +107,7 @@ interaction_matrix SH I N N N N N N N N N N I I I I I I N N N N N N N N N N N I sidechain_cutoffs default 3.0 4.0 # COO sidechain_cutoffs COO COO 2.5 3.5 -Sidechain_cutoffs COO SER 2.65 3.65 +sidechain_cutoffs COO SER 2.65 3.65 sidechain_cutoffs COO ARG 1.85 2.85 sidechain_cutoffs COO LYS 2.85 3.85 sidechain_cutoffs COO HIS 2.0 3.0 diff --git a/propka/protonate.py b/propka/protonate.py index 93306cb..49a1f8d 100644 --- a/propka/protonate.py +++ b/propka/protonate.py @@ -9,10 +9,16 @@ """ import logging import math +from typing import Iterable, TYPE_CHECKING + import propka.bonds import propka.atom +from propka.atom import Atom from propka.vector_algebra import rotate_vector_around_an_axis, Vector +if TYPE_CHECKING: + from propka.molecular_container import MolecularContainer + _LOGGER = logging.getLogger(__name__) @@ -58,7 +64,7 @@ def __init__(self, verbose=False): 'Br': 1.41, 'I': 1.61, 'S': 1.35} self.protonation_methods = {4: self.tetrahedral, 3: self.trigonal} - def protonate(self, molecules): + def protonate(self, molecules: "MolecularContainer"): """Protonate all atoms in the molecular container. Args: @@ -75,7 +81,7 @@ def protonate(self, molecules): self.protonate_atom(atom) @staticmethod - def remove_all_hydrogen_atoms(molecular_container): + def remove_all_hydrogen_atoms(molecular_container: "MolecularContainer"): """Remove all hydrogen atoms from molecule. Args: @@ -86,7 +92,7 @@ def remove_all_hydrogen_atoms(molecular_container): molecular_container.conformations[name] .get_non_hydrogen_atoms()) - def set_charge(self, atom): + def set_charge(self, atom: Atom): """Set charge for atom. Args: @@ -109,7 +115,7 @@ def set_charge(self, atom): atom.sybyl_type = atom.sybyl_type.replace('-', '') atom.charge_set = True - def protonate_atom(self, atom): + def protonate_atom(self, atom: Atom): """Protonate an atom. Args: @@ -126,7 +132,7 @@ def protonate_atom(self, atom): atom.is_protonated = True @staticmethod - def set_proton_names(heavy_atoms): + def set_proton_names(heavy_atoms: Iterable[Atom]): """Set names for protons. Args: @@ -139,7 +145,7 @@ def set_proton_names(heavy_atoms): bonded.name += str(i) i += 1 - def set_number_of_protons_to_add(self, atom): + def set_number_of_protons_to_add(self, atom: Atom): """Set the number of protons to add to this atom. Args: @@ -169,7 +175,7 @@ def set_number_of_protons_to_add(self, atom): _LOGGER.debug('-'*10) _LOGGER.debug(atom.number_of_protons_to_add) - def set_steric_number_and_lone_pairs(self, atom): + def set_steric_number_and_lone_pairs(self, atom: Atom): """Set steric number and lone pairs for atom. Args: @@ -207,8 +213,7 @@ def set_steric_number_and_lone_pairs(self, atom): atom.steric_number += 0 _LOGGER.debug('{0:>65s}: {1:>4.1f}'.format( 'Charge(-)', atom.charge)) - atom.steric_number -= atom.charge - atom.steric_number = math.floor(atom.steric_number/2.0) + atom.steric_number = math.floor((atom.steric_number - atom.charge) / 2) atom.number_of_lone_pairs = ( atom.steric_number - len(atom.bonded_atoms) - atom.number_of_protons_to_add @@ -220,7 +225,7 @@ def set_steric_number_and_lone_pairs(self, atom): 'Number of lone pairs', atom.number_of_lone_pairs)) atom.steric_num_lone_pairs_set = True - def add_protons(self, atom): + def add_protons(self, atom: Atom): """Add protons to atom. Args: @@ -236,7 +241,7 @@ def add_protons(self, atom): '(steric number: {0:d})'.format(atom.steric_number) ) - def trigonal(self, atom): + def trigonal(self, atom: Atom): """Add hydrogens in trigonal geometry. Args: @@ -269,15 +274,15 @@ def trigonal(self, atom): vec2 = Vector(atom1=atom.bonded_atoms[0], atom2=atom.bonded_atoms[0] .bonded_atoms[other_atom_indices[0]]) - axis = vec1**vec2 + axis = vec1.cross(vec2) # this is a trick to make sure that the order of atoms doesn't # influence the final postions of added protons if len(other_atom_indices) > 1: vec3 = Vector(atom1=atom.bonded_atoms[0], atom2=atom.bonded_atoms[0] .bonded_atoms[other_atom_indices[1]]) - axis2 = vec1**vec3 - if axis*axis2 > 0: + axis2 = vec1.cross(vec3) + if axis.dot(axis2) > 0: axis = axis+axis2 else: axis = axis-axis2 @@ -296,7 +301,7 @@ def trigonal(self, atom): new_a = self.set_bond_distance(new_a, atom.element) self.add_proton(atom, cvec+new_a) - def tetrahedral(self, atom): + def tetrahedral(self, atom: Atom): """Protonate atom in tetrahedral geometry. Args: @@ -338,13 +343,14 @@ def tetrahedral(self, atom): self.add_proton(atom, cvec+new_a) @staticmethod - def add_proton(atom, position): + def add_proton(atom: Atom, position: Vector): """Add a proton to an atom at a specific position. Args: atom: atom to protonate position: position for proton """ + assert atom.conformation_container is not None # Create the new proton new_h = propka.atom.Atom() new_h.set_property( @@ -368,7 +374,6 @@ def add_proton(atom, position): new_h.number_of_lone_pairs = 0 new_h.number_of_protons_to_add = 0 new_h.num_pi_elec_2_3_bonds = 0 - new_h.is_protonates = True atom.bonded_atoms.append(new_h) atom.number_of_protons_to_add -= 1 atom.conformation_container.add_atom(new_h) @@ -386,7 +391,7 @@ def add_proton(atom, position): i += 1 _LOGGER.debug('added %s %s %s', new_h, 'to', atom) - def set_bond_distance(self, bvec, element): + def set_bond_distance(self, bvec: Vector, element: str) -> Vector: """Set bond distance between atom and element. Args: diff --git a/propka/py.typed b/propka/py.typed new file mode 100644 index 0000000..e69de29 diff --git a/propka/run.py b/propka/run.py index ba8d53f..5e4bf65 100644 --- a/propka/run.py +++ b/propka/run.py @@ -12,10 +12,12 @@ """ import logging import sys +from typing import IO, Iterable, Optional from propka.lib import loadOptions from propka.input import read_parameter_file, read_molecule_file from propka.parameters import Parameters from propka.molecular_container import MolecularContainer +from propka.output import _PathArg _LOGGER = logging.getLogger(__name__) @@ -44,7 +46,9 @@ def main(optargs=None): my_molecule.write_pka() -def single(filename: str, optargs: tuple = (), stream=None, +def single(filename: _PathArg, + optargs: Iterable[str] = (), + stream: Optional[IO[str]] = None, write_pka: bool = True): """Run a single PROPKA calculation using ``filename`` as input. @@ -105,6 +109,7 @@ def single(filename: str, optargs: tuple = (), stream=None, .. versionchanged:: 3.4.0 Removed ability to write out PROPKA input files. """ + filename = str(filename) # Deal with input optarg options optargs = tuple(optargs) optargs += (filename,) diff --git a/propka/vector_algebra.py b/propka/vector_algebra.py index 1c65ab7..b216037 100644 --- a/propka/vector_algebra.py +++ b/propka/vector_algebra.py @@ -6,7 +6,8 @@ """ import logging import math -from typing import Optional, Protocol, Union +from typing import Optional, Protocol, overload +import warnings _LOGGER = logging.getLogger(__name__) @@ -69,20 +70,31 @@ def __sub__(self, other: _XYZ): self.y - other.y, self.z - other.z) - def __mul__(self, other: Union["Vector", "Matrix4x4", float]): + def dot(self, other: _XYZ) -> float: + return self.x * other.x + self.y * other.y + self.z * other.z + + @overload + def __mul__(self, other: "Vector") -> float: + ... + + @overload + def __mul__(self, other: "Matrix4x4") -> "Vector": + ... + + @overload + def __mul__(self, other: float) -> "Vector": + ... + + def __mul__(self, other): """Dot product, scalar and matrix multiplication.""" if isinstance(other, Vector): - return self.x * other.x + self.y * other.y + self.z * other.z - elif isinstance(other, Matrix4x4): - return Vector( - xi=other.a11*self.x + other.a12*self.y + other.a13*self.z - + other.a14*1.0, - yi=other.a21*self.x + other.a22*self.y + other.a23*self.z - + other.a24*1.0, - zi=other.a31*self.x + other.a32*self.y + other.a33*self.z - + other.a34*1.0 - ) - elif type(other) in [int, float]: + warnings.warn("Use Vector.dot() instead of operator.mul()", DeprecationWarning, stacklevel=2) + return self.dot(other) + if isinstance(other, Matrix4x4): + warnings.warn("Use M @ v (operator.matmul()) instead of M * v (operator.mul())", + DeprecationWarning, stacklevel=2) + return other @ self + if isinstance(other, (int, float)): return Vector(self.x * other, self.y * other, self.z * other) raise TypeError(f'{type(other)} not supported') @@ -90,6 +102,10 @@ def __rmul__(self, other): return self.__mul__(other) def __pow__(self, other: _XYZ): + warnings.warn("Use Vector.cross() instead of operator.pow()", DeprecationWarning, stacklevel=2) + return self.cross(other) + + def cross(self, other: _XYZ): """Cross product.""" return Vector(self.y * other.z - self.z * other.y, self.z * other.x - self.x * other.z, @@ -160,6 +176,17 @@ def __init__(self, self.a43 = a43i self.a44 = a44i + def __matmul__(self, v: _XYZ) -> Vector: + """Matrix vector multiplication with homogeneous coordinates. + + Assumes that the last row is (0, 0, 0, 1). + """ + return Vector( + self.a11 * v.x + self.a12 * v.y + self.a13 * v.z + self.a14, + self.a21 * v.x + self.a22 * v.y + self.a23 * v.z + self.a24, + self.a31 * v.x + self.a32 * v.y + self.a33 * v.z + self.a34, + ) + def angle(avec: Vector, bvec: Vector) -> float: """Get the angle between two vectors. @@ -170,7 +197,7 @@ def angle(avec: Vector, bvec: Vector) -> float: Returns: angle in radians """ - dot = avec * bvec + dot = avec.dot(bvec) return math.acos(dot / (avec.length() * bvec.length())) @@ -196,10 +223,10 @@ def signed_angle_around_axis(avec: Vector, bvec: Vector, axis: Vector) -> float: Returns: angle in radians """ - norma = avec**axis - normb = bvec**axis + norma = avec.cross(axis) + normb = bvec.cross(axis) ang = angle(norma, normb) - dot_ = bvec*(avec**axis) + dot_ = bvec.dot(avec.cross(axis)) if dot_ < 0: ang = -ang return ang @@ -223,21 +250,21 @@ def rotate_vector_around_an_axis(theta: float, axis: Vector, vec: Vector) -> Vec else: gamma = math.pi/2.0 rot_z = rotate_atoms_around_z_axis(gamma) - vec = rot_z * vec - axis = rot_z * axis + vec = rot_z @ vec + axis = rot_z @ axis beta = 0.0 if axis.x != 0: beta = -axis.x/abs(axis.x)*math.acos( axis.z/math.sqrt(axis.x*axis.x + axis.z*axis.z)) rot_y = rotate_atoms_around_y_axis(beta) - vec = rot_y * vec - axis = rot_y * axis + vec = rot_y @ vec + axis = rot_y @ axis rot_z = rotate_atoms_around_z_axis(theta) - vec = rot_z * vec + vec = rot_z @ vec rot_y = rotate_atoms_around_y_axis(-beta) - vec = rot_y * vec + vec = rot_y @ vec rot_z = rotate_atoms_around_z_axis(-gamma) - vec = rot_z * vec + vec = rot_z @ vec return vec diff --git a/propka/version.py b/propka/version.py index 9d1edf5..ef5cec2 100644 --- a/propka/version.py +++ b/propka/version.py @@ -7,6 +7,7 @@ TODO - this module unnecessarily confuses the code. Can we eliminate it? """ import logging +from propka.atom import Atom from propka.hydrogens import setup_bonding_and_protonation, setup_bonding from propka.hydrogens import setup_bonding_and_protonation_30_style from propka.energy import radial_volume_desolvation, calculate_pair_weight @@ -98,6 +99,10 @@ def setup_bonding(self, molecular_container): """Setup bonding using assigned model.""" return self.prepare_bonds(self.parameters, molecular_container) + def get_hydrogen_bond_parameters(self, atom1: Atom, atom2: Atom) -> tuple: + """Get hydrogen bond parameters for two atoms.""" + raise NotImplementedError("abstract method") + class VersionA(Version): """TODO - figure out what this is.""" diff --git a/tests/test_basic_regression.py b/tests/test_basic_regression.py index 8f27263..11b9d61 100644 --- a/tests/test_basic_regression.py +++ b/tests/test_basic_regression.py @@ -68,11 +68,11 @@ def run_propka(options, pdb_path, tmp_path): """ options += [str(pdb_path)] args = loadOptions(options) + cwd = Path.cwd() try: _LOGGER.warning( "Working in tmpdir {0:s} because of PROPKA file output; " "need to fix this.".format(str(tmp_path))) - cwd = Path.cwd() os.chdir(tmp_path) parameters = read_parameter_file(args.parameters, Parameters()) molecule = MolecularContainer(parameters, args) @@ -148,6 +148,7 @@ def compare_output(pdb, tmp_path, ref_path): def test_regression(pdb, options, tmp_path): """Basic regression test of PROPKA functionality.""" path_dict = get_test_dirs() + ref_path = None for ext in ["json", "dat"]: ref_path = path_dict["results"] / f"{pdb}.{ext}" diff --git a/tests/test_lib.py b/tests/test_lib.py new file mode 100644 index 0000000..8a9c66e --- /dev/null +++ b/tests/test_lib.py @@ -0,0 +1,27 @@ +import propka.lib as m + +import argparse +import pytest + + +def test_parse_res_string(): + assert m.parse_res_string("C:123") == ("C", 123, " ") + assert m.parse_res_string("C:123B") == ("C", 123, "B") + assert m.parse_res_string("ABC:123x") == ("ABC", 123, "x") + with pytest.raises(ValueError): + m.parse_res_string("C:B123") + with pytest.raises(ValueError): + m.parse_res_string("123B") + with pytest.raises(ValueError): + m.parse_res_string("C:123:B") + + +def test_parse_res_list(): + assert m.parse_res_list("C:123") == [("C", 123, " ")] + assert m.parse_res_list("ABC:123,D:4,F:56X") == [ + ("ABC", 123, " "), + ("D", 4, " "), + ("F", 56, "X"), + ] + with pytest.raises(argparse.ArgumentTypeError): + m.parse_res_list("C:B123") diff --git a/tests/test_vector_algebra.py b/tests/test_vector_algebra.py new file mode 100644 index 0000000..0cde8f0 --- /dev/null +++ b/tests/test_vector_algebra.py @@ -0,0 +1,163 @@ +import propka.vector_algebra as m + +import math +import pytest +from pytest import approx + +RADIANS_90 = math.pi / 2 + + +def assert_vector_equal(v1: m.Vector, v2: m.Vector): + assert isinstance(v1, m.Vector) + assert isinstance(v2, m.Vector) + assert v1.x == approx(v2.x) + assert v1.y == approx(v2.y) + assert v1.z == approx(v2.z) + + +def _matrix4x4_tolist(self: m.Matrix4x4) -> list: + return [ + self.a11, self.a12, self.a13, self.a14, self.a21, self.a22, self.a23, self.a24, + self.a31, self.a32, self.a33, self.a34, self.a41, self.a42, self.a43, self.a44 + ] + +def assert_matrix4x4_equal(m1: m.Matrix4x4, m2: m.Matrix4x4): + assert isinstance(m1, m.Matrix4x4) + assert isinstance(m2, m.Matrix4x4) + assert _matrix4x4_tolist(m1) == approx(_matrix4x4_tolist(m2)) + + +def test_Vector__init(): + v = m.Vector() + assert v.x == 0.0 + assert v.y == 0.0 + assert v.z == 0.0 + v = m.Vector(12, 34, 56) + assert v.x == 12 + assert v.y == 34 + assert v.z == 56 + v1 = m.Vector(atom1=v) + assert v1.x == 12 + assert v1.y == 34 + assert v1.z == 56 + v2 = m.Vector(5, 4, 3) + v3 = m.Vector(atom1=v2, atom2=v1) + assert v3.x == 7 + assert v3.y == 30 + assert v3.z == 53 + + +def test_Vector__plusminus(): + v1 = m.Vector(1, 2, 3) + v2 = m.Vector(4, 5, 6) + v3 = v1 + v2 + assert v3.x == 5 + assert v3.y == 7 + assert v3.z == 9 + v3 = v1 - v2 + assert v3.x == -3 + assert v3.y == -3 + assert v3.z == -3 + v4 = -v1 + assert v4.x == -1 + assert v4.y == -2 + assert v4.z == -3 + + +def test_Vector__mul__number(): + v1 = m.Vector(1, 2, 3) + assert_vector_equal(v1 * 2, m.Vector(2, 4, 6)) + + +def test_Vector__mul__Vector(): + v1 = m.Vector(1, 2, 3) + v2 = m.Vector(4, 5, 6) + with pytest.deprecated_call(): + assert v1 * v2 == 32 + assert v1.dot(v2) == 32 + with pytest.raises(TypeError): + v1 @ v2 # type: ignore + + +def test_Vector__mul__Matrix4x4(): + v1 = m.Vector(1, 2, 3) + assert_vector_equal(m.Matrix4x4() @ v1, m.Vector()) + m2 = m.Matrix4x4(0, 1, 0, 0, 20, 0, 0, 0, 0, 0, 300, 0, 0, 0, 0, 1) + with pytest.deprecated_call(): + assert_vector_equal(v1 * m2, m.Vector(2, 20, 900)) + with pytest.deprecated_call(): + assert_vector_equal(m2 * v1, m.Vector(2, 20, 900)) + assert_vector_equal(m2 @ v1, m.Vector(2, 20, 900)) + with pytest.raises(TypeError): + v1 @ m2 # type: ignore + + +def test_Vector__cross(): + v1 = m.Vector(1, 2, 3) + v2 = m.Vector(4, 5, 6) + with pytest.deprecated_call(): + assert_vector_equal(v1**v2, m.Vector(-3, 6, -3)) + assert_vector_equal(v1.cross(v2), m.Vector(-3, 6, -3)) + assert_vector_equal(v2.cross(v1), m.Vector(3, -6, 3)) + + +def test_Vector__length(): + v1 = m.Vector(1, 2, 3) + assert v1.length() == 14**0.5 + assert v1.sq_length() == 14 + + +def test_Vector__orthogonal(): + v1 = m.Vector(1, 2, 3) + assert v1.dot(v1.orthogonal()) == 0 + + +def test_Vector__rescale(): + v1 = m.Vector(1, 2, 3) + v2 = v1.rescale(4) + assert v2.length() == 4 + assert v2.x / v1.x == approx(4 / 14**0.5) + assert v2.y / v1.y == approx(4 / 14**0.5) + assert v2.z / v1.z == approx(4 / 14**0.5) + + +def test_angle(): + v1 = m.Vector(0, 0, 1) + v2 = m.Vector(0, 2, 0) + assert m.angle(v1, v2) == RADIANS_90 + + +def test_angle_degrees(): + v1 = m.Vector(0, 0, 3) + v2 = m.Vector(5, 0, 0) + assert m.angle_degrees(v1, v2) == 90 + + +def test_signed_angle_around_axis(): + v1 = m.Vector(0, 0, 3) + v2 = m.Vector(5, 0, 0) + v3 = m.Vector(0, 1, 0) + assert m.signed_angle_around_axis(v1, v2, v3) == -RADIANS_90 + v1 = m.Vector(0, 2, 3) + v2 = m.Vector(5, 4, 0) + assert m.signed_angle_around_axis(v1, v2, v3) == -RADIANS_90 + + +def test_rotate_vector_around_an_axis(): + v1 = m.Vector(0, 0, 3) + v2 = m.Vector(3, 0, 0) + v3 = m.Vector(0, -1, 0) + v4 = m.rotate_vector_around_an_axis(RADIANS_90, v3, v2) + assert_vector_equal(v4, v1) + + +def test_rotate_atoms_around_z_axis(): + m_rot = m.rotate_atoms_around_z_axis(-RADIANS_90) + assert_matrix4x4_equal(m_rot, + m.Matrix4x4(0, 1, 0, 0, -1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1)) + + +def test_rotate_atoms_around_y_axis(): + m_rot = m.rotate_atoms_around_y_axis(RADIANS_90) + assert_matrix4x4_equal(m_rot, + m.Matrix4x4(0, 0, 1, 0, 0, 1, 0, 0, -1, 0, 0, 0, 0, 0, 0, 1))