# # * This library is free software; you can redistribute it and/or # * modify it under the terms of the GNU Lesser General Public # * License as published by the Free Software Foundation; either # * version 2.1 of the License, or (at your option) any later version. # * # * This library is distributed in the hope that it will be useful, # * but WITHOUT ANY WARRANTY; without even the implied warranty of # * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU # * Lesser General Public License for more details. # #propka3.0, revision 182 2011-08-09 #------------------------------------------------------------------------------------------------------- #-- -- #-- PROPKA: A PROTEIN PKA PREDICTOR -- #-- -- #-- VERSION 3.0, 01/01/2011, COPENHAGEN -- #-- BY MATS H.M. OLSSON AND CHRESTEN R. SONDERGARD -- #-- -- #------------------------------------------------------------------------------------------------------- # # #------------------------------------------------------------------------------------------------------- # References: # # Very Fast Empirical Prediction and Rationalization of Protein pKa Values # Hui Li, Andrew D. Robertson and Jan H. Jensen # PROTEINS: Structure, Function, and Bioinformatics 61:704-721 (2005) # # Very Fast Prediction and Rationalization of pKa Values for Protein-Ligand Complexes # Delphine C. Bas, David M. Rogers and Jan H. Jensen # PROTEINS: Structure, Function, and Bioinformatics 73:765-783 (2008) # # PROPKA3: Consistent Treatment of Internal and Surface Residues in Empirical pKa predictions # Mats H.M. Olsson, Chresten R. Sondergard, Michal Rostkowski, and Jan H. Jensen # Journal of Chemical Theory and Computation, 7, 525-537 (2011) #------------------------------------------------------------------------------------------------------- import sys import string import math import copy from . import lib from .pdb import Atom pka_print = lib.pka_print class Residue: """ Residue class - contains atoms and properties of the residue """ def __init__(self, atoms, resName=None, resNumb=None, chainID=None, resInfo=None, options=None): """ Constructer of the residue object. """ self.label = None self.atoms = [] if chainID == None: self.chainID = atoms[0].chainID else: self.chainID = chainID if resNumb == None: self.resNumb = atoms[0].resNumb else: self.resNumb = resNumb if resName == None: self.resName = atoms[0].resName else: self.resName = resName self.resType = None # determins the interaction parameters, e.g. 'COO' self.Q = None # residue charge self.type = None # 'amino-acid' / 'ion' / 'ligand' / 'N-terminus' / 'C-terminus' self.buried = None self.location= None self.determinants = [[], [], []] self.Nmass = 0 self.Nlocl = 0 self.Emass = 0.0 self.Vmass = 0.0 self.Elocl = 0.0 self.Vlocl = 0.0 self.x = 0.00 self.y = 0.00 self.z = 0.00 self.pKa_mod = None self.pKa_pro = None self.pKas = [] # list with several pKa-objects: not used yet self.center_atoms = [] # list of atoms constituting the 'residue center', needed! self.default_key = None # the 'default' configuration if not all are there self.configurations = [] # keys to the atom configurations belonging to this residue self.coupled_residues = [] # setting residue information i.e. resType, type, Q, pKa_mod & pKa_pro from resInfo dictionary self.setResidueInformation(resInfo=resInfo) # setting residue label of the residue to create self.setResidueLabel() # setting up the atom labels used to calculate the 'residue center' residue_center_atom_labels = lib.residueCenterAtomList(self.resName) # adding/setting various atom-related properties such as atoms, configurations etc. for atom in atoms: # checking if this residue has 'HETATM' atoms; good chance it's a ligand then if atom.type == "hetatm" and self.type == None: self.type = "ligand" # setting the number of configurations for this residue for key in list(atom.configurations.keys()): if key not in self.configurations: self.configurations.append(key) #if atom.name[0] != 'H': if True: atom.set_residue(self) self.atoms.append(atom) # setting 'center atoms'; needs to be on self since it is reset when you switch configurations if len(residue_center_atom_labels) == 0: self.center_atoms.append(atom) elif atom.name in residue_center_atom_labels: self.center_atoms.append(atom) # set residue center, i.e. give 'x, y, z' values if len(self.center_atoms) > 0: self.setResidueCenter() def setResidueInformation(self, resInfo=None): """ set residue information based on resName - it is set here since it is a convenience thing """ # resType - determines interaction parameters if self.resName in resInfo['resType']: self.resType = resInfo['resType'][self.resName] # Q - group charge if self.resName in resInfo['Q']: self.Q = resInfo['Q'][self.resName] else: self.Q = 0.00 # type - 'amino-acid' / 'ion' / 'ligand' / 'N-terminus' / 'C-terminus' if self.resName in lib.residueList("standard"): self.type = "amino-acid" elif self.resName in resInfo['type']: self.type = resInfo['type'][self.resName] # pKa_mod - model or water pKa value if self.resName in resInfo['pKa']: self.pKa_mod = resInfo['pKa'][self.resName] self.pKa_pro = resInfo['pKa'][self.resName] else: self.pKa_mod = resInfo['pKa']['default'] self.pKa_pro = resInfo['pKa']['default'] def setConfiguration(self, key=None): """ set the 'current possition' to a specific 'configuration', or default if it doesn't exist """ self.cleanupPKA() if key in self.configurations: configuration = key else: configuration = self.default_key for atom in self.atoms: atom.setConfiguration(key=configuration) self.setResidueCenter() def cleanupPKA(self): """ Initializing/cleaning up residue! """ self.Nmass = 0 self.Nlocl = 0 self.Emass = 0.00 self.Elocl = 0.00 self.buried = 0.00 self.pKa_pro = self.pKa_mod self.determinants = [[], [], []] def setResidueCenter(self): """ sets the center of the residue based on center_atoms """ number_of_atoms = len(self.center_atoms) self.x = 0.00 self.y = 0.00 self.z = 0.00 for atom in self.center_atoms: self.x += atom.x self.y += atom.y self.z += atom.z if number_of_atoms > 0: self.x = self.x/number_of_atoms self.y = self.y/number_of_atoms self.z = self.z/number_of_atoms def getThirdAtomInAngle(self, atom=None): """ finds and returns the third atom in angular dependent interactions expecting one of ["HIS", "ARG", "AMD", "TRP"] """ if self.resName == "HIS": if atom.name == "HD1": return self.getAtom(name="ND1") elif atom.name == "HE2": return self.getAtom(name="NE2") elif self.resName == "ARG": if atom.name in ["HE"]: return self.getAtom(name="NE") elif atom.name in ["1HH1", "2HH1"]: return self.getAtom(name="NH1") elif atom.name in ["1HH2", "2HH2"]: return self.getAtom(name="NH2") elif self.resName == "ASN": return self.getAtom(name="ND2") elif self.resName == "GLN": return self.getAtom(name="NE2") elif self.resName == "TRP": return self.getAtom(name="NE1") def getAtom(self, name=None): """ finds and returns the specified atom in this residue. """ for atom in self.atoms: if atom.name == name: return atom break return None def checkOXT(self): """ Checks that OXT is present or creates it. """ O = self.getAtom(name="O") OXT = self.getAtom(name="OXT") # NMR Xplor over-write if O == None: O = self.getAtom(name="OT1") if O != None: O.setProperty(name="O") if OXT == None: OXT = self.getAtom(name="OT2") if OXT != None: OXT.setProperty(name="OXT") # continuing after 'NMR exception'; creating OXT if not found if OXT == None: # did not find OXT, creating it 'on the fly' CA = self.getAtom(name="CA") C = self.getAtom(name="C") if O == None or CA == None or C == None: pka_print("ERROR: cannot create OXT atom - missing CA, C, or O atoms; please correct pdbfile"); sys.exit(8) dX = -((CA.x-C.x) + (O.x-C.x)) dY = -((CA.y-C.y) + (O.y-C.y)) dZ = -((CA.z-C.z) + (O.z-C.z)) distance = math.sqrt( dX*dX + dY*dY + dZ*dZ ) x = C.x + 1.23*dX/distance y = C.y + 1.23*dY/distance z = C.z + 1.23*dZ/distance OXT = C.makeCopy(name="OXT", x=x, y=y, z=z) self.atoms.append(OXT) pka_print("creating %s atom" % (OXT.name)) return [O, OXT] def fillUnknownConfigurations(self, keys=None, options=None): """ Fills in the configurations that have not been read """ # getting default key as the first OK element in the sorted protein keys for key in keys: if key in self.configurations: self.default_key = key break # enforcing all residue keys on each atom for atom in self.atoms: for configuration in self.configurations: if key not in atom.configurations: atom.configurations[configuration] = atom.configurations[self.default_key] def checked(self, options=None): """ Checks that I understand all residues. """ excluded_resNames = lib.residueList("excluded") if self.resName in excluded_resNames: return False else: return True def checkResidue(self, options=None): """ Checks that I understand all residues. """ residue_list = lib.residueList("all") rename = {"LYP": "LYS", "CYX": "CYS", "HSD": "HIS", } # renaming residue if in rename dictionary if self.resName in rename: newName = rename[self.resName] if options.verbose == True: pka_print("Warning: renaming %s to %s" % (self.resName, newName)) self.resName = newName for atom in self.atoms: atom.resName = newName else: """OK - lets rock""" # after rename residue cases, check heavy atoms if self.resName in residue_list: # OK - lets rock self.checkAtoms(options=options) # Chresten's stuff elif self.type == "ion": outstr = "%s%4d - OK %s" % (self.resName, self.resNumb, self.type) pka_print(outstr) #elif self.resName in version.ions.keys(): # str = "%-3s%4d - %s with charge %+d" % (self.resName, # self.resNumb, # version.ions_long_names[self.resName], # version.ions[self.resName]) # pka_print(str) else: outstr = "%s%4d - unidentified residue" % (self.resName, self.resNumb) pka_print(outstr) def checkAtoms(self, options=None): """ Checks that all heavy atoms are there """ outstr = "%s%4d - " % (self.resName, self.resNumb) atom_list = lib.atomList(self.resName) OK = True for name in atom_list: FOUND = False for atom in self.atoms: if atom.name == name: FOUND = True if FOUND == False: outstr += " %s" % (name) OK = False if OK == True: self.checkConfigurations(verbose=False) outstr += " OK (%2d: %2d)" % (len(self.atoms), len(self.configurations)) if options.verbose == True: pka_print(outstr) else: outstr += " missing" pka_print(outstr) def checkConfigurations(self, verbose=False): """ checks that all atoms in this residue has the same number of configurations """ for atom in self.atoms: for key in self.configurations: if key not in atom.configurations: atom.configurations[key] = atom.configurations[self.default_key] def printLabel(self): """ prints the residue ID """ outstr = "%s%4d %s" % (self.resName, self.resNumb, self.chainID) pka_print(outstr) def __str__(self): return self.label def extractBackBoneAtoms(self): """ Returns the back-bone atoms """ N = None H = None C = None O = None for atom in self.atoms: if atom.name == "N": N = atom elif atom.name == "H": H = atom elif atom.name == "C": C = atom elif atom.name == "O": O = atom return N, H, C, O def makeDeterminantAtomList(self, resType=None, type=None, version=None): """ Extracting reference to determinant atom - test stage still """ if type == 'base' or type == 'acid': pair_type = type else: if resType in ["HIS", "LYS", "ARG", "N+ "]: pair_type = 'base' else: pair_type = 'acid' if self.resName not in version.atomInteractionList[pair_type]: pka_print("cannot find atomInteractionList for residue %s in residue.makeDeterminantAtomList()" % (self.resName)) sys.exit(9) # Searching for determinant atom atoms = [] for atom in self.atoms: if atom.name in version.atomInteractionList[pair_type][self.resName]: atoms.append(atom) return atoms def calculateTotalPKA(self): """ Calculates the total pKa values from the desolvation and determinants """ back_bone = 0.00 for determinant in self.determinants[0]: value = determinant.value back_bone += value side_chain = 0.00 for determinant in self.determinants[1]: value = determinant.value side_chain += value coulomb = 0.00 for determinant in self.determinants[2]: value = determinant.value coulomb += value self.pKa_pro = self.pKa_mod + self.Emass + self.Elocl + back_bone + side_chain + coulomb def calculateIntrinsicPKA(self): """ Calculates the intrinsic pKa values from the desolvation determinants, back-bone hydrogen bonds, and side-chain hydrogen bond to non-titratable residues """ back_bone = 0.00 for determinant in self.determinants[1]: value = determinant.value back_bone += value side_chain = 0.00 for determinant in self.determinants[0]: if determinant.label[0:3] not in ['ASP','GLU','LYS','ARG','HIS','CYS','TYR','C- ','N+ ']: value = determinant.value side_chain += value self.intrinsic_pKa = self.pKa_mod + self.Emass + self.Elocl + back_bone + side_chain return def calculateDesolvation(self, atoms, version=None, options=None): """ Calculates the desolvation contribution """ version.calculateDesolvation(self, atoms, options=options) def setChain(self, chainID): """ Set a chainID """ self.chainID = chainID for atom in self.atoms: atom.chainID = chainID def setResidueNumber(self, resNumb): """ Set the residue numbers to 'resNumb' """ self.resNumb = resNumb def setResidueLabel(self, label=None): """ Set the residue label to e.g. 'GLU 145 A' """ if label == None: self.label = lib.makeResidueLabel(self.resName, self.resNumb, self.chainID) else: self.label = label def shiftResidueNumber(self, shift): """ Shift the residue numbers with 'shift' """ self.resNumb = self.resNumb + shift def calculateFoldingEnergy(self, pH=None, reference=None, options=None): """ returning the electrostatic energy of this residue at pH 'pH' """ if pH == None: pH = options.pH if reference == None: reference = options.reference # calculating the ddG(neutral --> low-pH) contribution if self.resType not in ["COO", "HIS", "N+ ", "CYS", "TYR", "LYS", "ARG"]: ddG = 0.00 else: if reference == "low-pH": ddG_neutral = 0.00 else: if self.Q > 0.00: pKa_prime = self.pKa_pro coulomb_determinants = self.determinants[2] for determinant in coulomb_determinants: if determinant.value > 0.00: pKa_prime -= determinant.value ddG_neutral = -1.36*(pKa_prime - self.pKa_mod) else: ddG_neutral = 0.00 # calculating the ddG(low-pH --> pH) contribution # folded x = pH - self.pKa_pro y = 10**x Q_pro = math.log10(1+y) # unfolded x = pH - self.pKa_mod y = 10**x Q_mod = math.log10(1+y) ddG_low = -1.36*(Q_pro - Q_mod) ddG = ddG_neutral + ddG_low return ddG def getCharge(self, pH, state): """ returning the charge of this residue at pH 'pH' """ if state == "mod" or state == "unfolded": x = self.Q*(self.pKa_mod - pH) else: x = self.Q*(self.pKa_pro - pH) #x = pH - self.pKa_pro y = 10**x charge = self.Q*(y/(1.0+y)) #charge = math.log10(1+y) return charge def calculateTitrationCurve(self, grid=[0., 14., 0.10]): """ calculates the titration curve of this residue """ if grid == None: grid = [self.pKa_pro-2.5, self.pKa_pro+2.5, 0.10] state = "folded" titration_curve = [] pH = grid[0] stop = grid[1] + grid[2]/2.0 while pH < stop: Q = self.getCharge(pH, state) titration_curve.append([pH, Q]) pH += grid[2] return titration_curve def getSummaryString(self): """ Writing the summary string """ outstr = " %s%8.2lf%10.2lf" % (self.label, self.pKa_pro, self.pKa_mod) return outstr def mutateToAla(self): """ mutating residue to 'ALA' Note, if you mutate a ionizable residue to Ala it will remain in 'ionizable_residues list' and propka will try to calcualte the pKa of it !!! """ keep_atoms = lib.atomList("ALA") self.printLabel() self.resName = "ALA" self.setResidueLabel() self.printLabel() new_atoms = [] for atom in self.atoms: if atom.name in keep_atoms: new_atoms.append(atom) pka_print(self.atoms) pka_print(new_atoms) self.cleanupResidue self.pKa_mod = pKa_mod(self.resName) self.pKa_pro = self.pKa_mod self.atoms = new_atoms def replaceWithResidue(self, new_residue): """ replacing current residue with incoming residue """ self.cleanupResidue() self.resName = new_residue.resName self.pKa_mod = new_residue.pKa_mod self.setResidueLabel() exclude_atoms = ["N", "CA", "C", "O"] tmp_atoms = [] for atom in self.atoms: if atom.name in exclude_atoms: tmp_atoms.append(atom) self.atoms = tmp_atoms for new_atom in new_residue.atoms: if new_atom.name in exclude_atoms: """ do nothing """ else: self.atoms.append(new_atom) def getDeterminantString(self): """ Everything should be calculated, now, let's print the darn thing and be done! """ #if self.location == "SURFACE" or self.location == "BURIED ": BURIED_RATIO = True empty_determinant = "%s%4d%2s" % ("XXX", 0, "X") number_of_sidechain = len(self.determinants[0]) number_of_backbone = len(self.determinants[1]) number_of_coulomb = len(self.determinants[2]) number_of_determinants = number_of_sidechain + number_of_backbone + number_of_coulomb number_of_lines = max(1, number_of_sidechain, number_of_backbone, number_of_coulomb) outsting = "" #outsting += " number_of_sidechain = %d" % (number_of_sidechain) #outsting += " number_of_backbone = %d" % (number_of_backbone ) #outsting += " number_of_coulomb = %d" % (number_of_coulomb ) #outsting += " number_of_lines = %d" % (number_of_lines ) #print outsting for line_number in range(1, number_of_lines+1): outsting += "%s" % (self.label) if line_number == 1: outsting += " %6.2lf" % (self.pKa_pro) if len(self.coupled_residues)>0: outsting+='*' else: outsting+=' ' if BURIED_RATIO == True: if self.type == "BONDED": outsting += " BONDED " else: outsting += " %4d%2s " % (int(100.0*self.buried), "%") else: outsting += "%8s" % (self.type) outsting += " %6.2lf %4d" % (self.Emass, self.Nmass) outsting += " %6.2lf %4d" % (self.Elocl, self.Nlocl) else: outsting += "%40s" % (" ") # Side-chain determinants if line_number > number_of_sidechain : outsting += "%8.2lf %s" % (0.0, empty_determinant) else: determinant = self.determinants[0][line_number-1] outsting += "%8.2lf %s" % (determinant.value, determinant.label) # Back-bone determinants if line_number > number_of_backbone: outsting += "%8.2lf %s" % (0.0, empty_determinant) else: determinant = self.determinants[1][line_number-1] outsting += "%8.2lf %s" % (determinant.value, determinant.label) # Coulomb determinants if line_number > number_of_coulomb: outsting += "%8.2lf %s" % (0.0, empty_determinant) else: determinant = self.determinants[2][line_number-1] outsting += "%8.2lf %s" % (determinant.value, determinant.label) # adding end-of-line outsting += "\n" return outsting def printResult(self): """ Everything should be calculated, now, let's print the darn thing and be done! """ #if self.location == "SURFACE" or self.location == "BURIED ": BURIED_RATIO = True empty_determinant = "%s%4d%2s" % ("XXX", 0, "X") number_of_sidechain = len(self.determinants[0]) number_of_backbone = len(self.determinants[1]) number_of_coulomb = len(self.determinants[2]) number_of_determinants = number_of_sidechain + number_of_backbone + number_of_coulomb number_of_lines = max(1, number_of_sidechain, number_of_backbone, number_of_coulomb) outstr = "" outstr += " number_of_sidechain = %d" % (number_of_sidechain) outstr += " number_of_backbone = %d" % (number_of_backbone ) outstr += " number_of_coulomb = %d" % (number_of_coulomb ) outstr += " number_of_lines = %d" % (number_of_lines ) #print outstr if True: for line_number in range(1, number_of_lines+1): outstr = "%s%4d%2s" % (self.resName, self.resNumb, self.chainID) if line_number == 1: outstr += " %6.2lf" % (self.pKa_pro) if BURIED_RATIO == True: outstr += " %4d%2s " % (int(100.0*self.buried), "%") else: outstr += " %8s" % (self.type) outstr += " %6.2lf %4d" % (self.Emass, self.Nmass) outstr += " %6.2lf %4d" % (self.Elocl, self.Nlocl) else: outstr += "%40s" % (" ") # Side-chain determinant if line_number > number_of_sidechain : outstr += "%8.2lf %s" % (0.0, empty_determinant) else: determinant = self.determinants[0][line_number-1] outstr += "%8.2lf %s" % (determinant.value, determinant.label) # Back-bone determinant if line_number > number_of_backbone: outstr += "%8.2lf %s" % (0.0, empty_determinant) else: determinant = self.determinants[1][line_number-1] outstr += "%8.2lf %s" % (determinant.value, determinant.label) # Coulomb determinant if line_number > number_of_coulomb: outstr += "%8.2lf %s" % (0.0, empty_determinant) else: determinant = self.determinants[2][line_number-1] outstr += "%8.2lf %s" % (determinant.value, determinant.label) pka_print('%s' % (outstr)) else: outstr = "%s%4d%2s%4d%2d" % (self.resName, self.resNumb, self.chainID, number_of_lines, number_of_determinants) pka_print('%s' % (outstr)) pka_print('') def translate(self, translation): """ translate residue according to 'translation' """ for atom in self.atoms: atom.x += translation[0] atom.y += translation[1] atom.z += translation[2] for key in list(atom.configurations.keys()): for i in range(3): atom.configurations[key][i] += translation[i] def rotate(self, axis, theta, center=None): """ rotate residue theta radians around axis with center=center """ from .rotate import generalRotationMatrix translate = [0.00, 0.00, 0.00] number_of_atoms = 0 for atom in self.atoms: if atom.name in center or center == None: number_of_atoms += 1 translate[0] += atom.x/len(self.atoms) translate[1] += atom.y/len(self.atoms) translate[2] += atom.z/len(self.atoms) for atom in self.atoms: for i in range(3): translate[i] = translate[i]/number_of_atoms # translate to rotation center for atom in self.atoms: atom.x -= translate[0] atom.y -= translate[1] atom.z -= translate[2] # get rotation matrix rotation_matrix = generalRotationMatrix(axis, theta) # rotating new_position = [None, None, None] for atom in self.atoms: # rotate actual position old_position = [atom.x, atom.y, atom.z] for xyz in range(3): new_position[xyz] = translate[xyz] for i in range(3): new_position[xyz] += rotation_matrix[xyz][i]*old_position[i] # update position atom.x = new_position[0] atom.y = new_position[1] atom.z = new_position[2] # rotate configuration for key in list(atom.configurations.keys()): for xyz in range(3): new_position[xyz] = translate[xyz] for i in range(3): new_position[xyz] += rotation_matrix[xyz][i]*atom.configurations[key][i] for xyz in range(3): atom.configurations[key][xyz] = new_position[xyz] def makeCopy(self, chainID = None, resNumb = None, ): """ making a copy of this residue """ from .protein import getResidueParameters if chainID == None: chainID = self.chainID if resNumb == None: resNumb = self.resNumb newAtoms = [] for atom in self.atoms: newAtoms.append( atom.makeCopy(chainID=chainID, resNumb=resNumb) ) resInfo = getResidueParameters() newResidue = Residue(newAtoms, resInfo=resInfo) newResidue.resType = self.resType newResidue.Q = self.Q newResidue.type = self.type return newResidue