Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

De-lint PROPKA #33

Merged
merged 15 commits into from
May 23, 2020
Merged
Show file tree
Hide file tree
Changes from 8 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
538 changes: 281 additions & 257 deletions propka/atom.py

Large diffs are not rendered by default.

554 changes: 250 additions & 304 deletions propka/bonds.py

Large diffs are not rendered by default.

123 changes: 39 additions & 84 deletions propka/calculations.py
Original file line number Diff line number Diff line change
@@ -1,97 +1,52 @@

from __future__ import division
from __future__ import print_function

import math, propka.protonate, propka.bonds,copy, sys
"""PROPKA calculations."""
import math
import copy
import sys
import propka.protonate
import propka.bonds
from propka.lib import info, warning


#
# methods for setting up atom naming, bonding, and protonation
#

def setup_bonding_and_protonation(parameters, molecular_container):
"""Set up bonding and protonation for a molecule.

Args:
molecular_container: molecule container.
"""
# make bonds
my_bond_maker = setup_bonding(parameters, molecular_container)

# set up ligand atom names
set_ligand_atom_names(molecular_container)

# apply information on pi electrons
my_bond_maker.add_pi_electron_information(molecular_container)

# Protonate atoms
if molecular_container.options.protonate_all:
my_protonator = propka.protonate.Protonate(verbose=False)
my_protonator.protonate(molecular_container)


return



def setup_bonding(parameters, molecular_container):
# make bonds
my_bond_maker = propka.bonds.bondmaker()
my_bond_maker.find_bonds_for_molecules_using_boxes(molecular_container)
"""Set up bonding for a molecular container.

Args:
molecular_container: the molecular container in question
Returns:
BondMaker object
"""
my_bond_maker = propka.bonds.BondMaker()
my_bond_maker.find_bonds_for_molecules_using_boxes(molecular_container)
return my_bond_maker


def set_ligand_atom_names(molecular_container):
for name in molecular_container.conformation_names:
molecular_container.conformations[name].set_ligand_atom_names()
return



#
# The following methods imitates propka3.0 protonation!
#
def setup_bonding_and_protonation_30_style(parameters, molecular_container):
# Protonate atoms
protonate_30_style(molecular_container)
"""Set names for ligands in molecular container.

# make bonds
my_bond_maker = propka.bonds.bondmaker()
my_bond_maker.find_bonds_for_molecules_using_boxes(molecular_container)

return


def protonate_30_style(molecular_container):
Args:
molecular_container: molecular container for ligand names
"""
for name in molecular_container.conformation_names:
info('Now protonating', name)
# split atom into residues
curres = -1000000
residue = []
O=None
C=None
for atom in molecular_container.conformations[name].atoms:
if atom.resNumb != curres:
curres = atom.resNumb
if len(residue)>0:
#backbone
[O, C]= addBackBoneHydrogen(residue,O,C)
#arginine
if residue[0].resName == 'ARG':
addArgHydrogen(residue)
#histidine
if residue[0].resName == 'HIS':
addHisHydrogen(residue)
#tryptophan
if residue[0].resName == 'TRP':
addTrpHydrogen(residue)
#amides
if residue[0].resName in ['GLN','ASN']:
addAmdHydrogen(residue)


residue = []
if atom.type=='atom':
residue.append(atom)
molecular_container.conformations[name].set_ligand_atom_names()

return

def addArgHydrogen(residue):
"""
Expand Down Expand Up @@ -159,7 +114,7 @@ def addTrpHydrogen(residue):
elif atom.name == "CE2":
CE = atom
if CD == None or NE == None or CE == None:
str = "Did not find all atoms in %s%4d - in %s" % (self.resName, self.resNumb, "addTrpHydrogen()")
str = "Did not find all atoms in %s%4d - in %s" % (self.res_name, self.res_num, "addTrpHydrogen()")
info(str)
sys.exit(0)

Expand All @@ -176,15 +131,15 @@ def addAmdHydrogen(residue):
O = None
N = None
for atom in residue:
if (atom.resName == "GLN" and atom.name == "CD") or (atom.resName == "ASN" and atom.name == "CG"):
if (atom.res_name == "GLN" and atom.name == "CD") or (atom.res_name == "ASN" and atom.name == "CG"):
C = atom
elif (atom.resName == "GLN" and atom.name == "OE1") or (atom.resName == "ASN" and atom.name == "OD1"):
elif (atom.res_name == "GLN" and atom.name == "OE1") or (atom.res_name == "ASN" and atom.name == "OD1"):
O = atom
elif (atom.resName == "GLN" and atom.name == "NE2") or (atom.resName == "ASN" and atom.name == "ND2"):
elif (atom.res_name == "GLN" and atom.name == "NE2") or (atom.res_name == "ASN" and atom.name == "ND2"):
N = atom

if C == None or O == None or N == None:
str = "Did not find N, C and/or O in %s%4d - in %s" % (atom.resName, atom.resNumb, "addAmdHydrogen()")
str = "Did not find N, C and/or O in %s%4d - in %s" % (atom.res_name, atom.res_num, "addAmdHydrogen()")
info(str)
sys.exit(0)

Expand Down Expand Up @@ -222,7 +177,7 @@ def addBackBoneHydrogen(residue, O, C):
return [new_O,new_C]


if N.resName == "PRO":
if N.res_name == "PRO":
""" PRO doesn't have an H-atom; do nothing """
else:
H = protonateDirection([N, O, C])
Expand Down Expand Up @@ -311,11 +266,11 @@ def protonateSP2(list):
def make_new_H(atom, x,y,z):

new_H = propka.atom.Atom()
new_H.setProperty(numb = None,
new_H.set_property(numb = None,
name = 'H%s'%atom.name[1:],
resName = atom.resName,
chainID = atom.chainID,
resNumb = atom.resNumb,
res_name = atom.res_name,
chain_id = atom.chain_id,
res_num = atom.res_num,
x = x,
y = y,
z = z,
Expand All @@ -329,7 +284,7 @@ def make_new_H(atom, x,y,z):
new_H.steric_number = 0
new_H.number_of_lone_pairs = 0
new_H.number_of_protons_to_add = 0
new_H.number_of_pi_electrons_in_double_and_triple_bonds = 0
new_H.num_pi_elec_2_3_bonds = 0

atom.bonded_atoms.append(new_H)
atom.conformation_container.add_atom(new_H)
Expand Down Expand Up @@ -357,7 +312,7 @@ def radial_volume_desolvation(parameters, group):

for atom in all_atoms:
# ignore atoms in the same residue
if atom.resNumb == group.atom.resNumb and atom.chainID == group.atom.chainID:
if atom.res_num == group.atom.res_num and atom.chain_id == group.atom.chain_id:
continue

sq_dist = squared_distance(group, atom)
Expand Down Expand Up @@ -409,15 +364,15 @@ def contactDesolvation(parameters, group):
'N+': 4.5}

all_atoms = group.atom.conformation_container.get_non_hydrogen_atoms()
if residue.resName in version.desolvationRadii:
local_cutoff = version.desolvationRadii[residue.resName]
if residue.res_name in version.desolvationRadii:
local_cutoff = version.desolvationRadii[residue.res_name]
else:
local_cutoff = 0.00
residue.Nmass = 0
residue.Nlocl = 0

for atom in all_atoms:
if atom.resNumb != group.atom.resNumb or atom.chainID != group.atom.chainID:
if atom.res_num != group.atom.res_num or atom.chain_id != group.atom.chain_id:
dX = atom.x - residue.x
dY = atom.y - residue.y
dZ = atom.z - residue.z
Expand Down
14 changes: 7 additions & 7 deletions propka/conformation_container.py
Original file line number Diff line number Diff line change
Expand Up @@ -158,7 +158,7 @@ def init_group(self, group):
titrate_only = self.molecular_container.options.titrate_only
if titrate_only is not None:
at = group.atom
if not (at.chainID, at.resNumb, at.icode) in titrate_only:
if not (at.chain_id, at.res_num, at.icode) in titrate_only:
group.titratable = False
if group.residue_type == 'CYS':
group.exclude_cys_from_results = True
Expand Down Expand Up @@ -376,7 +376,7 @@ def get_heavy_ligand_atoms(self):
return [atom for atom in self.atoms if atom.type=='hetatm' and atom.element != 'H']

def get_chain(self,chain):
return [atom for atom in self.atoms if atom.chainID != chain]
return [atom for atom in self.atoms if atom.chain_id != chain]


def add_atom(self, atom):
Expand All @@ -388,13 +388,13 @@ def add_atom(self, atom):
atom.molecular_container = self.molecular_container

# store chain id for bookkeeping
if not atom.chainID in self.chains:
self.chains.append(atom.chainID)
if not atom.chain_id in self.chains:
self.chains.append(atom.chain_id)

return

def copy_atom(self, atom):
new_atom = atom.makeCopy()
new_atom = atom.make_copy()
self.atoms.append(new_atom)
new_atom.conformation_container = self

Expand Down Expand Up @@ -443,8 +443,8 @@ def sort_atoms(self):
return

def sort_atoms_key(self, atom):
key = ord(atom.chainID)*1e7
key += atom.resNumb*1000
key = ord(atom.chain_id)*1e7
key += atom.res_num*1000
if len(atom.name) > len(atom.element):
key += ord(atom.name[len(atom.element)])
#info(atom,ord(atom.name[len(atom.element)]), '|%s||%s|'%(atom.name,atom.element))
Expand Down
2 changes: 1 addition & 1 deletion propka/determinants.py
Original file line number Diff line number Diff line change
Expand Up @@ -54,7 +54,7 @@ def addDeterminants(group1, group2, distance, version):
def addSidechainDeterminants(group1, group2, version=None):
"""
adding side-chain determinants/perturbations
Note, resNumb1 > resNumb2
Note, res_num1 > res_num2
"""

hbond_interaction = version.hydrogen_bond_interaction(group1, group2)
Expand Down
26 changes: 13 additions & 13 deletions propka/group.py
Original file line number Diff line number Diff line change
Expand Up @@ -111,22 +111,22 @@ def __init__(self, atom):
self.common_charge_centre = False


self.residue_type = self.atom.resName
self.residue_type = self.atom.res_name
if self.atom.terminal:
self.residue_type = self.atom.terminal

if self.atom.type=='atom':
self.label = '%-3s%4d%2s'%(self.residue_type, atom.resNumb, atom.chainID)
elif self.atom.resName in ['DA ','DC ','DG ','DT ']:
self.label = '%-3s%4d%2s'%(self.residue_type, atom.res_num, atom.chain_id)
elif self.atom.res_name in ['DA ','DC ','DG ','DT ']:
self.label = '%1s%1s%1s%4d%2s'%(self.residue_type[1],
atom.element,
atom.name.replace('\'','')[-1],
atom.resNumb,
atom.chainID)
atom.res_num,
atom.chain_id)

# self.label = '%1s%1s%1s%4d%2s'%(self.residue_type[1], atom.element,atom.name[-1], atom.resNumb, atom.chainID)
# self.label = '%1s%1s%1s%4d%2s'%(self.residue_type[1], atom.element,atom.name[-1], atom.res_num, atom.chain_id)
else:
self.label = '%-3s%4s%2s'%(self.residue_type, atom.name, atom.chainID)
self.label = '%-3s%4s%2s'%(self.residue_type, atom.name, atom.chain_id)


# container for squared distances
Expand Down Expand Up @@ -253,7 +253,7 @@ def __eq__(self, other):
return self.label==other.label
else:
# For heterogene atoms we also need to check the residue number
return self.label==other.label and self.atom.resNumb == other.atom.resNumb
return self.label==other.label and self.atom.res_num == other.atom.res_num

def __hash__(self):
""" Needed together with __eq__ - otherwise we can't make sets of groups """
Expand Down Expand Up @@ -361,7 +361,7 @@ def setup(self):
if not self.model_pka_set:
self.model_pka = self.parameters.model_pkas[self.residue_type]
# check if we should apply a custom model pka
key = '%s-%s'%(self.atom.resName.strip(), self.atom.name.strip())
key = '%s-%s'%(self.atom.res_name.strip(), self.atom.name.strip())
if key in self.parameters.custom_model_pkas.keys():
self.model_pka = self.parameters.custom_model_pkas[key]

Expand Down Expand Up @@ -1122,7 +1122,7 @@ class Ion_group(Group):
def __init__(self, atom):
Group.__init__(self,atom)
self.type = 'ION'
self.residue_type = atom.resName.strip()
self.residue_type = atom.res_name.strip()
info('Found ion group:', atom)
return

Expand Down Expand Up @@ -1201,15 +1201,15 @@ def is_protein_group(parameters,atom):
### Backbone
if atom.type == 'atom' and atom.name == 'N':
# ignore proline backbone nitrogens
if atom.resName != 'PRO':
if atom.res_name != 'PRO':
return BBN_group(atom)
if atom.type == 'atom' and atom.name == 'C':
# ignore C- carboxyl
if atom.count_bonded_elements('O') == 1:
return BBC_group(atom)

### Filters for side chains based on PDB protein atom names
key = '%s-%s'%(atom.resName, atom.name)
key = '%s-%s'%(atom.res_name, atom.name)

if key in parameters.protein_group_mapping.keys():
return eval('%s_group(atom)'%parameters.protein_group_mapping[key])
Expand Down Expand Up @@ -1342,7 +1342,7 @@ def is_ligand_group_by_marvin_pkas(parameters, atom):

def is_ion_group(parameters, atom):

if atom.resName.strip() in parameters.ions.keys():
if atom.res_name.strip() in parameters.ions.keys():
return Ion_group(atom)

return None
12 changes: 6 additions & 6 deletions propka/iterative.py
Original file line number Diff line number Diff line change
Expand Up @@ -124,9 +124,9 @@ def addIterativeIonPair(object1, object2, interaction, version):
Q2 = object2.Q
comp1 = object1.pKa_old + annihilation[0] + Q1*coulomb_value
comp2 = object2.pKa_old + annihilation[1] + Q2*coulomb_value
if object1.resName not in version.parameters.exclude_sidechain_interactions:
if object1.res_name not in version.parameters.exclude_sidechain_interactions:
comp1 += Q1*hbond_value
if object2.resName not in version.parameters.exclude_sidechain_interactions:
if object2.res_name not in version.parameters.exclude_sidechain_interactions:
comp2 += Q2*hbond_value

if Q1 == -1.0 and comp1 < comp2:
Expand Down Expand Up @@ -155,12 +155,12 @@ def addIterativeIonPair(object1, object2, interaction, version):
# Side-chain
if hbond_value > 0.005:
# residue1
if object1.resName not in version.parameters.exclude_sidechain_interactions:
if object1.res_name not in version.parameters.exclude_sidechain_interactions:
interaction = [object2, Q1*hbond_value]
annihilation[0] += -Q1*hbond_value
object1.determinants['sidechain'].append(interaction)
# residue2
if object2.resName not in version.parameters.exclude_sidechain_interactions:
if object2.res_name not in version.parameters.exclude_sidechain_interactions:
interaction = [object1, Q2*hbond_value]
annihilation[1] += -Q2*hbond_value
object2.determinants['sidechain'].append(interaction)
Expand Down Expand Up @@ -309,7 +309,7 @@ def __init__(self, group):

self.label = group.label
self.atom = group.atom
self.resName = group.residue_type
self.res_name = group.residue_type
self.Q = group.charge
self.pKa_old = None
self.pKa_new = None
Expand Down Expand Up @@ -357,7 +357,7 @@ def __eq__(self, other):
return self.label==other.label
else:
# For heterogene atoms we also need to check the residue number
return self.label==other.label and self.atom.resNumb == other.atom.resNumb
return self.label==other.label and self.atom.res_num == other.atom.res_num

def __hash__(self):
""" Needed together with __eq__ - otherwise we can't make sets of groups """
Expand Down
Loading