1# 2# * This library is free software; you can redistribute it and/or 3# * modify it under the terms of the GNU Lesser General Public 4# * License as published by the Free Software Foundation; either 5# * version 2.1 of the License, or (at your option) any later version. 6# * 7# * This library is distributed in the hope that it will be useful, 8# * but WITHOUT ANY WARRANTY; without even the implied warranty of 9# * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 10# * Lesser General Public License for more details. 11# 12 13#propka3.0, revision 182 2011-08-09 14#------------------------------------------------------------------------------------------------------- 15#-- -- 16#-- PROPKA: A PROTEIN PKA PREDICTOR -- 17#-- -- 18#-- VERSION 3.0, 01/01/2011, COPENHAGEN -- 19#-- BY MATS H.M. OLSSON AND CHRESTEN R. SONDERGARD -- 20#-- -- 21#------------------------------------------------------------------------------------------------------- 22# 23# 24#------------------------------------------------------------------------------------------------------- 25# References: 26# 27# Very Fast Empirical Prediction and Rationalization of Protein pKa Values 28# Hui Li, Andrew D. Robertson and Jan H. Jensen 29# PROTEINS: Structure, Function, and Bioinformatics 61:704-721 (2005) 30# 31# Very Fast Prediction and Rationalization of pKa Values for Protein-Ligand Complexes 32# Delphine C. Bas, David M. Rogers and Jan H. Jensen 33# PROTEINS: Structure, Function, and Bioinformatics 73:765-783 (2008) 34# 35# PROPKA3: Consistent Treatment of Internal and Surface Residues in Empirical pKa predictions 36# Mats H.M. Olsson, Chresten R. Sondergard, Michal Rostkowski, and Jan H. Jensen 37# Journal of Chemical Theory and Computation, 7, 525-537 (2011) 38#------------------------------------------------------------------------------------------------------- 39import math, time 40 41import iterative 42import lib 43from lib import pka_print 44#import debug 45import calculator as calculate 46from determinant import Determinant 47 48 49def setDeterminants(propka_residues, version=None, options=None): 50 """ 51 adding side-chain and coulomb determinants/perturbations to all residues - note, backbone determinants are set separately 52 """ 53 #debug.printResidues(propka_residues) 54 iterative_interactions = [] 55 # --- NonIterative section ---# 56 for residue1 in propka_residues: 57 for residue2 in propka_residues: 58 59 if residue1 == residue2: 60 break 61 62 distance = calculate.InterResidueDistance(residue1, residue2) 63 64 if distance < version.sidechain_cutoff or distance < version.coulomb_cutoff[1]: 65 try: 66 do_pair, iterative_interaction = version.interaction[residue1.resType][residue2.resType] 67 except KeyError as e: 68 do_pair = False 69 70 if do_pair == True: 71 if iterative_interaction == True: 72 iterative.addtoDeterminantList(residue1, residue2, distance, iterative_interactions, version=version) 73 #print "%s - %s I" % (residue1.label, residue2.label) 74 else: 75 addDeterminants(residue1, residue2, distance, version=version) 76 #print "%s - %s" % (residue1.label, residue2.label) 77 else: 78 """ False - don't do this at home folks """ 79 80 # --- Iterative section ---# 81 #debug.printIterativeDeterminants(iterative_interactions) 82 iterative.addDeterminants(iterative_interactions, version, options=options) 83 84 85def addDeterminants(residue1, residue2, distance, version=None): 86 """ 87 adding determinants/perturbations, distance(R1, R2) < coulomb_cutoff always 88 """ 89 90 # side-chain determinant 91 if distance < version.sidechain_cutoff: 92 # Currently we don't want any hydrogen bonds to ligands 93 if 'ligand' not in [residue1.type, residue2.type]: 94 addSidechainDeterminants(residue1, residue2, version) 95 96 do_coulomb = version.checkCoulombPair(residue1, residue2, distance) 97 98 # Coulomb determinant 99 if do_coulomb == True: 100 addCoulombDeterminants(residue1, residue2, distance, version) 101 102 103def addSidechainDeterminants(residue1, residue2, version=None): 104 """ 105 adding side-chain determinants/perturbations 106 Note, resNumb1 > resNumb2 107 """ 108 distance = 999.0 109 closest_atom1 = None 110 closest_atom2 = None 111 atoms1 = residue1.makeDeterminantAtomList(residue2.resName, version=version) 112 atoms2 = residue2.makeDeterminantAtomList(residue1.resName, version=version) 113 for atom1 in atoms1: 114 for atom2 in atoms2: 115 # select the smallest inter-atom distance 116 current_distance = calculate.InterAtomDistance(atom1, atom2) 117 if current_distance < distance: 118 closest_atom1 = atom1 119 closest_atom2 = atom2 120 distance = current_distance 121 122 dpka_max, cutoff = version.SideChainParameters[residue1.resType][residue2.resType] 123 if distance < cutoff[1]: 124 if residue2.resType in version.angularDependentSideChainInteractions: 125 atom3 = residue2.getThirdAtomInAngle(closest_atom2) 126 distance, f_angle, nada = calculate.AngleFactorX(closest_atom1, closest_atom2, atom3) 127 elif residue1.resType in version.angularDependentSideChainInteractions: 128 atom3 = residue1.getThirdAtomInAngle(closest_atom1) 129 distance, f_angle, nada = calculate.AngleFactorX(closest_atom2, closest_atom1, atom3) 130 else: 131 # i.e. no angular dependence 132 f_angle = 1.0 133 134 weight = version.calculatePairWeight(residue1.Nmass, residue2.Nmass) 135 exception, value = version.checkExceptions(residue1, residue2) 136 #exception = False # circumventing exception 137 if exception == True: 138 """ do nothing, value should have been assigned """ 139 #pka_print(" exception for %s %s %6.2lf" % (residue1.label, residue2.label, value)) 140 else: 141 value = version.calculateSideChainEnergy(distance, dpka_max, cutoff, weight, f_angle) 142 if residue1.Q == residue2.Q: 143 # acid pair or base pair 144 if residue1.pKa_mod < residue2.pKa_mod: 145 newDeterminant1 = Determinant(residue2.label, -value) 146 newDeterminant2 = Determinant(residue1.label, value) 147 else: 148 newDeterminant1 = Determinant(residue2.label, value) 149 newDeterminant2 = Determinant(residue1.label, -value) 150 else: 151 newDeterminant1 = Determinant(residue2.label, value*residue1.Q) 152 newDeterminant2 = Determinant(residue1.label, value*residue2.Q) 153 if residue1.resName not in version.exclude_sidechain_interactions: 154 residue1.determinants[0].append(newDeterminant1) 155 if residue2.resName not in version.exclude_sidechain_interactions: 156 residue2.determinants[0].append(newDeterminant2) 157 158 159def addCoulombDeterminants(residue1, residue2, distance, version): 160 """ 161 adding NonIterative Coulomb determinants/perturbations 162 """ 163 weight = version.calculatePairWeight(residue1.Nmass, residue2.Nmass) 164 value = version.calculateCoulombEnergy(distance, weight) 165 Q1 = residue1.Q 166 Q2 = residue2.Q 167 168 # assigning the Coulombic interaction 169 if Q1 < 0.0 and Q2 < 0.0: 170 """ both are acids """ 171 addCoulombAcidPair(residue1, residue2, value) 172 elif Q1 > 0.0 and Q2 > 0.0: 173 """ both are bases """ 174 addCoulombBasePair(residue1, residue2, value) 175 else: 176 """ one of each """ 177 addCoulombIonPair(residue1, residue2, value) 178 179 180def addCoulombAcidPair(object1, object2, value): 181 """ 182 Adding the Coulomb interaction (an acid pair): 183 the higher pKa is raised 184 """ 185 label1 = object1.label 186 label2 = object2.label 187 if object1.pKa_mod > object2.pKa_mod: 188 newDeterminant = Determinant(label2, value) 189 object1.determinants[2].append(newDeterminant) 190 else: 191 newDeterminant = Determinant(label1, value) 192 object2.determinants[2].append(newDeterminant) 193 194 195def addCoulombBasePair(object1, object2, value): 196 """ 197 Adding the Coulomb interaction (a base pair): 198 the lower pKa is lowered 199 """ 200 label1 = object1.label 201 label2 = object2.label 202 if object1.pKa_mod < object2.pKa_mod: 203 newDeterminant = Determinant(label2, -value) 204 object1.determinants[2].append(newDeterminant) 205 else: 206 newDeterminant = Determinant(label1, -value) 207 object2.determinants[2].append(newDeterminant) 208 209 210def addCoulombIonPair(object1, object2, value): 211 """ 212 Adding the Coulomb interaction (an acid-base pair): 213 the pKa of the acid is lowered & the pKa of the base is raised 214 """ 215 label1 = object1.label 216 label2 = object2.label 217 218 # residue1 219 Q1 = object1.Q 220 newDeterminant = Determinant(label2, Q1*value) 221 object1.determinants[2].append(newDeterminant) 222 223 # residue2 224 Q2 = object2.Q 225 newDeterminant = Determinant(label1, Q2*value) 226 object2.determinants[2].append(newDeterminant) 227 228 229 230 231def setIonDeterminants(protein, version=None): 232 """ 233 adding ion determinants/perturbations 234 """ 235 ionizable_residues = lib.residueList("propka1") 236 for residue in protein.propka_residues: 237 if residue.resName in ionizable_residues: 238 for ion in protein.residue_dictionary["ION"]: 239 distance = calculate.InterResidueDistance(residue, ion) 240 if distance < version.coulomb_cutoff[1]: 241 label = "%s%4d%2s" % (ion.resName, ion.resNumb, ion.chainID) 242 weight = version.calculatePairWeight(residue.Nmass, ion.Nmass) 243 # the pKa of both acids and bases are shifted up by negative ions (and vice versa) 244 value = (-ion.Q) * version.calculateCoulombEnergy(distance, weight) 245 newDeterminant = Determinant(label, value) 246 residue.determinants[2].append(newDeterminant) 247 248 249def setBackBoneDeterminants(backbone_interactions, version=None): 250 """ 251 adding back-bone determinants/perturbations 252 Angle: atom1 -- atom2-atom3 253 backbone_interactions = [[acids, NH], [bases, CO]] 254 changing the code with minimum effect of method calls 255 """ 256 setBackBoneAcidDeterminants(backbone_interactions[0], version=version) 257 setBackBoneBaseDeterminants(backbone_interactions[1], version=version) 258 259 260def setBackBoneAcidDeterminants(data_clump, version=None): 261 """ 262 adding back-bone determinants/perturbations for acids: 263 Angle: atom1 -- atom2-atom3, i.e. COO -- H-N 264 data_clump = [acids, NH] 265 """ 266 residues, interactions = data_clump 267 for residue in residues: 268 if residue.location != "BONDED": 269 dpKa_max, cutoff = version.BackBoneParameters[residue.resType] 270 for interaction in interactions: 271 atom2 = interaction[1] 272 atom3 = interaction[0] 273 atoms = residue.makeDeterminantAtomList("back-bone", version=version) 274 shortest_distance = 999. 275 for atom in atoms: 276 distance = calculate.InterAtomDistance(atom, atom2) 277 if distance < shortest_distance: 278 shortest_distance = distance 279 atom1 = atom 280 distance, f_angle, nada = calculate.AngleFactorX(atom1, atom2, atom3) 281 if distance < cutoff[1] and f_angle > 0.001: 282 label = "%s%4d%2s" % (atom2.resName, atom2.resNumb, atom2.chainID) 283 value = residue.Q * calculate.HydrogenBondEnergy(distance, dpKa_max, cutoff, f_angle) 284 newDeterminant = Determinant(label, value) 285 residue.determinants[1].append(newDeterminant) 286 287 288def setBackBoneBaseDeterminants(data_clump, version=None): 289 """ 290 adding back-bone determinants/perturbations for bases: 291 Angle: atom1 -- atom2-atom3, i.e. C=O -- H-N(HIS) 292 data_clump = [bases, CO] 293 """ 294 residues, interactions = data_clump 295 for residue in residues: 296 if residue.location != "BONDED": 297 dpKa_max, cutoff = version.BackBoneParameters[residue.resType] 298 for interaction in interactions: 299 distance = 999. 300 atom1 = interaction[1] 301 atoms = residue.makeDeterminantAtomList("back-bone", version=version) 302 for atom in atoms: 303 current_distance = calculate.InterAtomDistance(atom1, atom) 304 if current_distance < distance: 305 atom2 = atom 306 distance = current_distance 307 if distance < cutoff[1]: 308 if residue.resType in version.angularDependentSideChainInteractions: 309 atom3 = residue.getThirdAtomInAngle(atom2) 310 distance, f_angle, nada = calculate.AngleFactorX(atom1, atom2, atom3) 311 else: 312 f_angle = 1.0 313 if f_angle > 0.001: 314 # add determinant 315 label = "%s%4d%2s" % (atom2.resName, atom2.resNumb, atom2.chainID) 316 value = residue.Q * calculate.HydrogenBondEnergy(distance, dpKa_max, cutoff, f_angle) 317 newDeterminant = Determinant(label, value) 318 residue.determinants[1].append(newDeterminant) 319 320