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