diff --git a/INSTALL b/INSTALL new file mode 100644 index 0000000..3865f52 --- /dev/null +++ b/INSTALL @@ -0,0 +1,19 @@ +## Installation + +Clone repository or unpack the tar ball and install with +[setuptools](http://pythonhosted.org/setuptools/index.html) (note: if +you don't have setuptools installed you will need an internet +connection so that the installation procedure can sownload the +required files): + + cd propka-3.1 + python setup.py install --user + +This will install the `propka31` script in your executable directory, +as configured for setuptools, for instance `~/.local/bin`. You can +change the bin directory with the `--install-scripts` option. For +example, in order to install in my `bin` directory in my home +directory: + + python setup.py install --user --install-scripts ~/bin + diff --git a/MANIFEST.in b/MANIFEST.in new file mode 100644 index 0000000..eb21e09 --- /dev/null +++ b/MANIFEST.in @@ -0,0 +1,4 @@ +include README.md INSTALL +include ez_setup.py setup.py +include propka.cfg + diff --git a/README.md b/README.md index 2864741..6e36338 100644 --- a/README.md +++ b/README.md @@ -1,45 +1,69 @@ # PROPKA 3.1 -PROPKA predicts the pKa values of ionizable groups in -proteins (version 3.0) and -protein-ligand complexes (version 3.1) -based in the 3D structure. +PROPKA predicts the pKa values of ionizable groups in proteins +(version 3.0) and protein-ligand complexes (version 3.1) +based on the 3D structure. For proteins without ligands both version should produce the same result. The method is described in the following papers, which you should cite in publications: -* Søndergaard, Chresten R., Mats HM Olsson, Michal Rostkowski, and Jan H. Jensen. "Improved Treatment of Ligands and Coupling Effects in Empirical Calculation and Rationalization of pKa Values." Journal of Chemical Theory and Computation 7, no. 7 (2011): 2284-2295. +* Sondergaard, Chresten R., Mats HM Olsson, Michal Rostkowski, and Jan H. Jensen. "Improved Treatment of Ligands and Coupling Effects in Empirical Calculation and Rationalization of pKa Values." Journal of Chemical Theory and Computation 7, no. 7 (2011): 2284-2295. -* Olsson, Mats HM, Chresten R. Søndergaard, Michal Rostkowski, and Jan H. Jensen. "PROPKA3: consistent treatment of internal and surface residues in empirical pKa predictions." Journal of Chemical Theory and Computation 7, no. 2 (2011): 525-537. +* Olsson, Mats HM, Chresten R. Sondergaard, Michal Rostkowski, and Jan H. Jensen. "PROPKA3: consistent treatment of internal and surface residues in empirical pKa predictions." Journal of Chemical Theory and Computation 7, no. 2 (2011): 525-537. See [propka.ki.ku.dk](http://propka.ki.ku.dk/) for the PROPKA web server, using the [tutorial](http://propka.ki.ku.dk/~luca/wiki/index.php/PROPKA_3.1_Tutorial). +## Modifications + +This release of PROPKA 3.1 was modified by Oliver Beckstein + from the released version. + +* Included patches from + https://github.com/schrodinger/propka-3.1/tree/python27-compat to + make it compatible with Python 2.7 + +* Packaged for installation with setuptools. + + ## Installation -No installation needed. Just clone and run. +Clone repository or unpack the tar ball and install with +[setuptools](http://pythonhosted.org/setuptools/index.html) (note: if +you don't have setuptools installed you will need an internet +connection so that the installation procedure can download the +required files): + + cd propka-3.1 + python setup.py install --user + +This will install the `propka31` script in your executable directory, +as configured for setuptools, for instance `~/.local/bin`. You can +change the bin directory with the `--install-scripts` option. For +example, in order to install in my `bin` directory in my home +directory: + + python setup.py install --user --install-scripts ~/bin + ## Requirements -* Python 3.1 or higher +* Python 2.7 or higher or Python 3.1 or higher ## Getting started 1. Clone the code from GitHub -2. Run 'propka.py' with a .pdb file (see Examples) +2. `python setup.py install --user` +2. Run `propka31` with a .pdb file (see Examples) ## Examples Calculate using pdb file - ./propka.py 1hpx.pdb + propka31 1hpx.pdb -If for some reason your setup with python3.1+ is -not located in '/usr/bin/python3', run the script - - python3.2 propka.py 1hpx.pdb ## Testing (for developers) @@ -49,9 +73,9 @@ Please run `Tests/runtest.py/` after changes before pushing commits. Please cite these references in publications: -* Søndergaard, Chresten R., Mats HM Olsson, Michal Rostkowski, and Jan H. Jensen. "Improved Treatment of Ligands and Coupling Effects in Empirical Calculation and Rationalization of pKa Values." Journal of Chemical Theory and Computation 7, no. 7 (2011): 2284-2295. +* Sondergaard, Chresten R., Mats HM Olsson, Michal Rostkowski, and Jan H. Jensen. "Improved Treatment of Ligands and Coupling Effects in Empirical Calculation and Rationalization of pKa Values." Journal of Chemical Theory and Computation 7, no. 7 (2011): 2284-2295. -* Olsson, Mats HM, Chresten R. Søndergaard, Michal Rostkowski, and Jan H. Jensen. "PROPKA3: consistent treatment of internal and surface residues in empirical pKa predictions." Journal of Chemical Theory and Computation 7, no. 2 (2011): 525-537. +* Olsson, Mats HM, Chresten R. Sondergaard, Michal Rostkowski, and Jan H. Jensen. "PROPKA3: consistent treatment of internal and surface residues in empirical pKa predictions." Journal of Chemical Theory and Computation 7, no. 2 (2011): 525-537. diff --git a/Tests/runtest.py b/Tests/runtest.py index 1ce29ab..e73517a 100755 --- a/Tests/runtest.py +++ b/Tests/runtest.py @@ -1,9 +1,13 @@ -#!/usr/bin/python3 +#!/usr/bin/env python """ Run test for test pdbs """ +from __future__ import division +from __future__ import print_function + from subprocess import call import os, re +import sys pdbs = ['1FTJ-Chain-A', '1HPX', @@ -14,7 +18,7 @@ for pdb in pdbs: print('RUNNING '+pdb) # Run pka calculation - call(['../propka.py','pdb/'+pdb+'.pdb'], stdout = open(pdb+'.out', 'w+')) + call([sys.executable, '../scripts/propka31.py','pdb/'+pdb+'.pdb'], stdout = open(pdb+'.out', 'w+')) # Test pka predictiona result = open('results/'+pdb+'.dat','r') diff --git a/ez_setup.py b/ez_setup.py new file mode 100644 index 0000000..2535472 --- /dev/null +++ b/ez_setup.py @@ -0,0 +1,258 @@ +#!python +"""Bootstrap setuptools installation + +If you want to use setuptools in your package's setup.py, just include this +file in the same directory with it, and add this to the top of your setup.py:: + + from ez_setup import use_setuptools + use_setuptools() + +If you want to require a specific version of setuptools, set a download +mirror, or use an alternate download directory, you can do so by supplying +the appropriate options to ``use_setuptools()``. + +This file can also be run as a script to install or upgrade setuptools. +""" +import os +import shutil +import sys +import tempfile +import tarfile +import optparse +import subprocess + +from distutils import log + +try: + from site import USER_SITE +except ImportError: + USER_SITE = None + +DEFAULT_VERSION = "0.9.6" +DEFAULT_URL = "https://pypi.python.org/packages/source/s/setuptools/" + +def _python_cmd(*args): + args = (sys.executable,) + args + return subprocess.call(args) == 0 + +def _install(tarball, install_args=()): + # extracting the tarball + tmpdir = tempfile.mkdtemp() + log.warn('Extracting in %s', tmpdir) + old_wd = os.getcwd() + try: + os.chdir(tmpdir) + tar = tarfile.open(tarball) + _extractall(tar) + tar.close() + + # going in the directory + subdir = os.path.join(tmpdir, os.listdir(tmpdir)[0]) + os.chdir(subdir) + log.warn('Now working in %s', subdir) + + # installing + log.warn('Installing Setuptools') + if not _python_cmd('setup.py', 'install', *install_args): + log.warn('Something went wrong during the installation.') + log.warn('See the error message above.') + # exitcode will be 2 + return 2 + finally: + os.chdir(old_wd) + shutil.rmtree(tmpdir) + + +def _build_egg(egg, tarball, to_dir): + # extracting the tarball + tmpdir = tempfile.mkdtemp() + log.warn('Extracting in %s', tmpdir) + old_wd = os.getcwd() + try: + os.chdir(tmpdir) + tar = tarfile.open(tarball) + _extractall(tar) + tar.close() + + # going in the directory + subdir = os.path.join(tmpdir, os.listdir(tmpdir)[0]) + os.chdir(subdir) + log.warn('Now working in %s', subdir) + + # building an egg + log.warn('Building a Setuptools egg in %s', to_dir) + _python_cmd('setup.py', '-q', 'bdist_egg', '--dist-dir', to_dir) + + finally: + os.chdir(old_wd) + shutil.rmtree(tmpdir) + # returning the result + log.warn(egg) + if not os.path.exists(egg): + raise IOError('Could not build the egg.') + + +def _do_download(version, download_base, to_dir, download_delay): + egg = os.path.join(to_dir, 'setuptools-%s-py%d.%d.egg' + % (version, sys.version_info[0], sys.version_info[1])) + if not os.path.exists(egg): + tarball = download_setuptools(version, download_base, + to_dir, download_delay) + _build_egg(egg, tarball, to_dir) + sys.path.insert(0, egg) + import setuptools + setuptools.bootstrap_install_from = egg + + +def use_setuptools(version=DEFAULT_VERSION, download_base=DEFAULT_URL, + to_dir=os.curdir, download_delay=15): + # making sure we use the absolute path + to_dir = os.path.abspath(to_dir) + was_imported = 'pkg_resources' in sys.modules or \ + 'setuptools' in sys.modules + try: + import pkg_resources + except ImportError: + return _do_download(version, download_base, to_dir, download_delay) + try: + pkg_resources.require("setuptools>=" + version) + return + except pkg_resources.VersionConflict: + e = sys.exc_info()[1] + if was_imported: + sys.stderr.write( + "The required version of setuptools (>=%s) is not available,\n" + "and can't be installed while this script is running. Please\n" + "install a more recent version first, using\n" + "'easy_install -U setuptools'." + "\n\n(Currently using %r)\n" % (version, e.args[0])) + sys.exit(2) + else: + del pkg_resources, sys.modules['pkg_resources'] # reload ok + return _do_download(version, download_base, to_dir, + download_delay) + except pkg_resources.DistributionNotFound: + return _do_download(version, download_base, to_dir, + download_delay) + + +def download_setuptools(version=DEFAULT_VERSION, download_base=DEFAULT_URL, + to_dir=os.curdir, delay=15): + """Download setuptools from a specified location and return its filename + + `version` should be a valid setuptools version number that is available + as an egg for download under the `download_base` URL (which should end + with a '/'). `to_dir` is the directory where the egg will be downloaded. + `delay` is the number of seconds to pause before an actual download + attempt. + """ + # making sure we use the absolute path + to_dir = os.path.abspath(to_dir) + try: + from urllib.request import urlopen + except ImportError: + from urllib2 import urlopen + tgz_name = "setuptools-%s.tar.gz" % version + url = download_base + tgz_name + saveto = os.path.join(to_dir, tgz_name) + src = dst = None + if not os.path.exists(saveto): # Avoid repeated downloads + try: + log.warn("Downloading %s", url) + src = urlopen(url) + # Read/write all in one block, so we don't create a corrupt file + # if the download is interrupted. + data = src.read() + dst = open(saveto, "wb") + dst.write(data) + finally: + if src: + src.close() + if dst: + dst.close() + return os.path.realpath(saveto) + + +def _extractall(self, path=".", members=None): + """Extract all members from the archive to the current working + directory and set owner, modification time and permissions on + directories afterwards. `path' specifies a different directory + to extract to. `members' is optional and must be a subset of the + list returned by getmembers(). + """ + import copy + import operator + from tarfile import ExtractError + directories = [] + + if members is None: + members = self + + for tarinfo in members: + if tarinfo.isdir(): + # Extract directories with a safe mode. + directories.append(tarinfo) + tarinfo = copy.copy(tarinfo) + tarinfo.mode = 448 # decimal for oct 0700 + self.extract(tarinfo, path) + + # Reverse sort directories. + if sys.version_info < (2, 4): + def sorter(dir1, dir2): + return cmp(dir1.name, dir2.name) + directories.sort(sorter) + directories.reverse() + else: + directories.sort(key=operator.attrgetter('name'), reverse=True) + + # Set correct owner, mtime and filemode on directories. + for tarinfo in directories: + dirpath = os.path.join(path, tarinfo.name) + try: + self.chown(tarinfo, dirpath) + self.utime(tarinfo, dirpath) + self.chmod(tarinfo, dirpath) + except ExtractError: + e = sys.exc_info()[1] + if self.errorlevel > 1: + raise + else: + self._dbg(1, "tarfile: %s" % e) + + +def _build_install_args(options): + """ + Build the arguments to 'python setup.py install' on the setuptools package + """ + install_args = [] + if options.user_install: + if sys.version_info < (2, 6): + log.warn("--user requires Python 2.6 or later") + raise SystemExit(1) + install_args.append('--user') + return install_args + +def _parse_args(): + """ + Parse the command line for options + """ + parser = optparse.OptionParser() + parser.add_option( + '--user', dest='user_install', action='store_true', default=False, + help='install in user site package (requires Python 2.6 or later)') + parser.add_option( + '--download-base', dest='download_base', metavar="URL", + default=DEFAULT_URL, + help='alternative URL from where to download the setuptools package') + options, args = parser.parse_args() + # positional arguments are ignored + return options + +def main(version=DEFAULT_VERSION): + """Install or upgrade setuptools and EasyInstall""" + options = _parse_args() + tarball = download_setuptools(download_base=options.download_base) + return _install(tarball, _build_install_args(options)) + +if __name__ == '__main__': + sys.exit(main()) diff --git a/propka.py b/propka.py deleted file mode 100755 index a789017..0000000 --- a/propka.py +++ /dev/null @@ -1,23 +0,0 @@ -#!/usr/bin/python3 - -import string,re,sys,os,math -import Source.lib, Source.molecular_container - -def main(): - """ - Reads in structure files, calculates pKa values, and prints pKa files - """ - # loading options, flaggs and arguments - options, pdbfiles = Source.lib.loadOptions() - - for pdbfile in pdbfiles: - my_molecule = Source.molecular_container.Molecular_container(pdbfile, options) - my_molecule.calculate_pka() - my_molecule.write_pka() - - -if __name__ == '__main__': - #import cProfile - #cProfile.run('main()',sort=1) - main() - diff --git a/Source/__init__.py b/propka/__init__.py similarity index 65% rename from Source/__init__.py rename to propka/__init__.py index e286beb..923c164 100644 --- a/Source/__init__.py +++ b/propka/__init__.py @@ -1,3 +1,5 @@ -# +# propka 3.1 +# (Module renamed from original Source to propka to facilitate setuptools installation) + __all__ = ['coupled_residues', 'lib', 'parameters', 'residue', 'bonds', 'debug', 'ligand', 'pdb', 'calculator', 'determinants', 'mutate', 'protein', 'chain', 'iterative', 'output', 'protonate'] diff --git a/Source/atom.py b/propka/atom.py similarity index 94% rename from Source/atom.py rename to propka/atom.py index c964930..b9abec4 100644 --- a/Source/atom.py +++ b/propka/atom.py @@ -1,4 +1,8 @@ -import string, Source.lib, Source.group + +from __future__ import division +from __future__ import print_function + +import string, propka.lib, propka.group class Atom: @@ -8,17 +12,17 @@ class Atom: def __init__(self, line=None, verbose=False): - self.set_properties(line) + self.set_properties(line) - self.residue_label = "%-3s%4d%2s" % (self.name,self.resNumb, self.chainID) + self.residue_label = "%-3s%4d%2s" % (self.name,self.resNumb, self.chainID) self.groups_extracted = 0 self.group = None self.group_type = None self.number_of_bonded_elements = {} - + self.cysteine_bridge = False - + self.bonded_atoms = [] self.residue = None self.conformation_container = None @@ -39,7 +43,7 @@ class Atom: self.sybyl_type = '' self.sybyl_assigned = False self.marvin_pka = False - + return @@ -122,13 +126,13 @@ class Atom: return True return False - + def setProperty(self, - numb = None, - name = None, - resName = None, + numb = None, + name = None, + resName = None, chainID = None, resNumb = None, x = None, @@ -157,7 +161,7 @@ class Atom: """ making a copy of this atom """ - + newAtom = Atom() newAtom.type = self.type newAtom.numb = self.numb @@ -191,11 +195,11 @@ class Atom: if self.group.titratable: model_pka = '%6.2f'%self.group.model_pka - + str = "%-6s%5d %s %s%2s%4d%12.3lf%8.3lf%8.3lf%6s%6s \n" % (self.type.upper(), self.numb, - Source.lib.makeTidyAtomLabel(self.name, + propka.lib.makeTidyAtomLabel(self.name, self.element), self.resName, self.chainID, @@ -213,7 +217,7 @@ class Atom: def make_conect_line(self): res = 'CONECT%5d'%self.numb - + # extract and sort numbers of bonded residues bonded = [] for atom in self.bonded_atoms: @@ -226,9 +230,9 @@ class Atom: res += '\n' return res - + def get_input_parameters(self): - """ Method for getting the input parameters stored in the + """ Method for getting the input parameters stored in the occupancy and b-factor fields in input files""" # Set the group type @@ -258,7 +262,7 @@ class Atom: # try to initialise the group try: - exec('self.group = Source.group.%s_group(self)'%self.occ) + exec('self.group = propka.group.%s_group(self)'%self.occ) except: raise Exception('%s in input_file is not recognized as a group'%self.occ) @@ -270,7 +274,7 @@ class Atom: # set occ and beta to standard values self.occ = '1.00' self.beta = '0.00' - + return @@ -282,7 +286,7 @@ class Atom: # making string str = "%-6s%5d %s %s%2s%4d%12.3lf%8.3lf%8.3lf%6s%6s\n" % (self.type.upper(), self.numb, - Source.lib.makeTidyAtomLabel(self.name, + propka.lib.makeTidyAtomLabel(self.name, self.element), self.resName, self.chainID, @@ -292,17 +296,17 @@ class Atom: self.z, self.occ, self.beta) - + return str def make_mol2_line(self,id): -#1 S1 3.6147 2.0531 1.4795 S.3 1 noname -0.1785 +#1 S1 3.6147 2.0531 1.4795 S.3 1 noname -0.1785 # making string str = "%-4d %-4s %10.4f %10.4f %10.4f %6s %6d %10s %10.4f\n" % (id, - Source.lib.makeTidyAtomLabel(self.name, + propka.lib.makeTidyAtomLabel(self.name, self.element), self.x, self.y, @@ -311,7 +315,7 @@ class Atom: self.resNumb, self.resName, 0.0)#self.charge) - + return str @@ -346,7 +350,7 @@ class Atom: # making string str = "ATOM " str += "%6d" % (numb) - str += " %s" % (Source.lib.makeTidyAtomLabel(name,self.element)) + str += " %s" % (propka.lib.makeTidyAtomLabel(name,self.element)) str += " %s" % (resName) str += "%2s" % (chainID) str += "%4d" % (resNumb) @@ -363,25 +367,25 @@ class Atom: """ Returns a 'tidier' atom label for printing the new pdbfile """ - return Source.lib.makeTidyAtomLabel(self.name,self.element) + return propka.lib.makeTidyAtomLabel(self.name,self.element) def __str__(self): return '%5d-%4s %5d-%3s (%1s) [%8.3f %8.3f %8.3f] %s' %(self.numb, self.name, self.resNumb, self.resName, self.chainID, self.x, self.y, self.z,self.element) - + # def get_element(self): # """ try to extract element if not already done""" # if self.element == '': # self.element = self.name.strip(string.digits)[0] # return self.element - + def set_residue(self, residue): """ Makes a references to the parent residue""" if self.residue == None: self.residue = residue - + diff --git a/Source/bonds.py b/propka/bonds.py similarity index 94% rename from Source/bonds.py rename to propka/bonds.py index bd3e72f..c0f7e71 100644 --- a/Source/bonds.py +++ b/propka/bonds.py @@ -1,14 +1,19 @@ -import pickle,sys,os,math,Source.calculations +from __future__ import division +from __future__ import print_function + +import pickle,sys,os,math,propka.calculations + +import pkg_resources class bondmaker: def __init__(self): - + # predefined bonding distances self.distances = { 'S-S':2.5, 'F-F':1.7} - + self.distances_squared = {} for k in self.distances.keys(): self.distances_squared[k]=self.distances[k]*self.distances[k] @@ -18,12 +23,11 @@ class bondmaker: self.H_dist_squared = self.H_dist * self.H_dist self.default_dist_squared = self.default_dist * self.default_dist - - self.max_sq_distance = max(list(self.distances_squared.values())+[self.default_dist_squared]) - + + self.max_sq_distance = max(list(self.distances_squared.values())+[self.default_dist_squared]) + # protein bonding data - path = os.path.split(__file__)[0] - self.data_file_name = os.path.join(path,'protein_bonds.dat') + self.data_file_name = pkg_resources.resource_filename(__name__, 'protein_bonds.dat') data = open(self.data_file_name,'rb') self.protein_bonds = pickle.load(data) @@ -104,13 +108,13 @@ class bondmaker: def find_bonds_for_protein(self, protein): - """ Finds bonds proteins based on the way atoms + """ Finds bonds proteins based on the way atoms normally bond in proteins""" print('++++ Side chains ++++') # side chains for chain in protein.chains: - for residue in chain.residues: + for residue in chain.residues: if residue.resName.replace(' ','') not in ['N+','C-']: self.find_bonds_for_side_chain(residue.atoms) @@ -145,7 +149,7 @@ class bondmaker: if atom1.name == 'SG': for atom2 in cys2.atoms: if atom2.name == 'SG': - if Source.calculations.squared_distance(atom1,atom2) < self.SS_dist_squared: + if propka.calculations.squared_distance(atom1,atom2) < self.SS_dist_squared: self.make_bond(atom1, atom2) @@ -174,7 +178,7 @@ class bondmaker: if atom1.name == 'C': for atom2 in residue2.atoms: if atom2.name == 'N': - if Source.calculations.squared_distance(atom1,atom2) < self.default_dist_squared: + if propka.calculations.squared_distance(atom1,atom2) < self.default_dist_squared: self.make_bond(atom1, atom2) return @@ -209,7 +213,7 @@ class bondmaker: for atom2 in atoms: if atom2.name in self.protein_bonds[atom1.resName][atom1.name]: self.make_bond(atom1,atom2) - + return @@ -270,7 +274,7 @@ class bondmaker: if atoms[i] in atoms[j].bonded_atoms: continue - + if self.check_distance(atoms[i], atoms[j]): self.make_bond(atoms[i],atoms[j]) # di-sulphide bonds @@ -283,14 +287,14 @@ class bondmaker: def check_distance(self, atom1, atom2): - sq_dist = Source.calculations.squared_distance(atom1, atom2) - + sq_dist = propka.calculations.squared_distance(atom1, atom2) + if sq_dist > self.max_sq_distance: return False - + key = '%s-%s'%(atom1.element,atom2.element) h_count = key.count('H') - + if sq_dist < self.H_dist_squared and h_count==1: return True if sq_dist < self.default_dist_squared and h_count==0: @@ -299,8 +303,8 @@ class bondmaker: if sq_dist < self.distances_squared[key]: return True - - + + return False def find_bonds_for_molecules_using_boxes(self, molecules): @@ -320,7 +324,7 @@ class bondmaker: def find_bonds_for_atoms_using_boxes(self, atoms): """ Finds all bonds for a list of atoms""" - + box_size = 2.5 # find min and max coordinates @@ -348,11 +352,13 @@ class bondmaker: #print('x range: [%6.2f;%6.2f] %6.2f'%(xmin,xmax,xlen)) #print('y range: [%6.2f;%6.2f] %6.2f'%(ymin,ymax,ylen)) #print('z range: [%6.2f;%6.2f] %6.2f'%(zmin,zmax,zlen)) - + # how many boxes do we need in each dimension? - self.no_box_x = max(1, math.ceil(xlen/box_size)) - self.no_box_y = max(1, math.ceil(ylen/box_size)) - self.no_box_z = max(1, math.ceil(zlen/box_size)) + # NOTE: math.ceil() returns an int in python3 and a float in python2, + # so we need to cast it to int for range() to work. + self.no_box_x = max(1, int(math.ceil(xlen/box_size))) + self.no_box_y = max(1, int(math.ceil(ylen/box_size))) + self.no_box_z = max(1, int(math.ceil(zlen/box_size))) #print('No. box x: %6.2f'%self.no_box_x) #print('No. box y: %6.2f'%self.no_box_y) @@ -372,12 +378,12 @@ class bondmaker: z = math.floor((atom.z-zmin)/box_size) self.put_atom_in_box(x,y,z,atom) - # assign bonds + # assign bonds keys = self.boxes.keys() for key in keys: self.find_bonds_for_atoms(self.boxes[key]) - + return @@ -386,7 +392,7 @@ class bondmaker: # one side of the x,y,z box in each dimension for bx in [x,x+1]: for by in [y,y+1]: - for bz in [z,z+1]: + for bz in [z,z+1]: key = self.box_key(bx,by,bz) if key in self.boxes.keys(): self.boxes[key].append(atom) @@ -397,13 +403,13 @@ class bondmaker: def box_key(self, x, y, z): return '%d-%d-%d'%(x,y,z) - + def has_bond(self, atom1, atom2): if atom1 in atom2.bonded_atoms or atom2 in atom1.bonded_atoms: return True return False - + def make_bond(self, atom1, atom2): """ Makes a bond between atom1 and atom2 """ @@ -413,13 +419,13 @@ class bondmaker: #print('making bond for',atom1,atom2) if not atom1 in atom2.bonded_atoms: atom2.bonded_atoms.append(atom1) - + if not atom2 in atom1.bonded_atoms: - atom1.bonded_atoms.append(atom2) + atom1.bonded_atoms.append(atom2) return - + def generate_protein_bond_dictionary(self, atoms): for atom in atoms: @@ -433,25 +439,25 @@ class bondmaker: not name_j in self.backbone_atoms: if not name_i in self.terminal_oxygen_names and\ not name_j in self.terminal_oxygen_names: - + if not resi_i in list(self.protein_bonds.keys()): self.protein_bonds[resi_i] = {} if not name_i in self.protein_bonds[resi_i]: self.protein_bonds[resi_i][name_i] = [] - + if not name_j in self.protein_bonds[resi_i][name_i]: self.protein_bonds[resi_i][name_i].append(name_j) - + if not resi_j in list(self.protein_bonds.keys()): self.protein_bonds[resi_j] = {} if not name_j in self.protein_bonds[resi_j]: self.protein_bonds[resi_j][name_j] = [] - + if not name_i in self.protein_bonds[resi_j][name_j]: self.protein_bonds[resi_j][name_j].append(name_i) - + return @@ -467,10 +473,10 @@ if __name__ == '__main__': if not os.path.isfile(filename): print('Error: Could not find \"%s\"'%filename) sys.exit(1) - + pdblist = pdb.readPDB(filename) my_protein = protein.Protein(pdblist,'test.pdb') - + for chain in my_protein.chains: for residue in chain.residues: residue.atoms = [atom for atom in residue.atoms if atom.element != 'H'] diff --git a/Source/calculations.py b/propka/calculations.py similarity index 94% rename from Source/calculations.py rename to propka/calculations.py index ed9f6aa..5f71f8e 100644 --- a/Source/calculations.py +++ b/propka/calculations.py @@ -1,6 +1,8 @@ +from __future__ import division +from __future__ import print_function -import math, Source.protonate, Source.bonds,copy, sys +import math, propka.protonate, propka.bonds,copy, sys # @@ -17,9 +19,9 @@ def setup_bonding_and_protonation(parameters, molecular_container): # apply information on pi electrons my_bond_maker.add_pi_electron_information(molecular_container) - # Protonate atoms + # Protonate atoms if molecular_container.options.protonate_all: - my_protonator = Source.protonate.Protonate(verbose=False) + my_protonator = propka.protonate.Protonate(verbose=False) my_protonator.protonate(molecular_container) @@ -29,7 +31,7 @@ def setup_bonding_and_protonation(parameters, molecular_container): def setup_bonding(parameters, molecular_container): # make bonds - my_bond_maker = Source.bonds.bondmaker() + my_bond_maker = propka.bonds.bondmaker() my_bond_maker.find_bonds_for_molecules_using_boxes(molecular_container) return my_bond_maker @@ -46,11 +48,11 @@ def set_ligand_atom_names(molecular_container): # The following methods imitates propka3.0 protonation! # def setup_bonding_and_protonation_30_style(parameters, molecular_container): - # Protonate atoms + # Protonate atoms protonate_30_style(molecular_container) # make bonds - my_bond_maker = Source.bonds.bondmaker() + my_bond_maker = propka.bonds.bondmaker() my_bond_maker.find_bonds_for_molecules_using_boxes(molecular_container) return @@ -87,7 +89,7 @@ def protonate_30_style(molecular_container): residue = [] if atom.type=='atom': residue.append(atom) - + return def addArgHydrogen(residue): @@ -117,7 +119,7 @@ def addArgHydrogen(residue): H4.name = "HN3" H5 = protonateDirection([NH2, NE, H1]) H5.name = "HN4" - + return [H1,H2,H3,H4,H5] def addHisHydrogen(residue): @@ -179,7 +181,7 @@ def addAmdHydrogen(residue): O = atom elif (atom.resName == "GLN" and atom.name == "NE2") or (atom.resName == "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()") print(str) @@ -189,20 +191,20 @@ def addAmdHydrogen(residue): H1.name = "HN1" H2 = protonateAverageDirection([N, C, O]) H2.name = "HN2" - + return def addBackBoneHydrogen(residue, O, C): """ - Adds hydrogen backbone atoms to residues according to the old way; dR is wrong for the N-terminus + Adds hydrogen backbone atoms to residues according to the old way; dR is wrong for the N-terminus (i.e. first residue) but it doesn't affect anything at the moment. Could be improved, but works for now. """ - + new_C = None new_O = None N = None - + for atom in residue: if atom.name == "N": @@ -211,20 +213,20 @@ def addBackBoneHydrogen(residue, O, C): new_C = atom if atom.name == "O": new_O = atom - - - + + + if None in [C,O,N]: return [new_O,new_C] - + if N.resName == "PRO": """ PRO doesn't have an H-atom; do nothing """ else: H = protonateDirection([N, O, C]) H.name = "H" - + return [new_O,new_C] @@ -247,7 +249,7 @@ def protonateDirection(list): y = X1.y + dY/length z = X1.z + dZ/length - + H = make_new_H(X1,x,y,z) H.name = "H" @@ -277,7 +279,7 @@ def protonateAverageDirection(list): H = make_new_H(X1,x,y,z) H.name = "H" - + return H @@ -290,7 +292,7 @@ def protonateSP2(list): X1 = list[0] X2 = list[1] X3 = list[2] - + dX = (X1.x + X3.x)*0.5 - X2.x dY = (X1.y + X3.y)*0.5 - X2.y dZ = (X1.z + X3.z)*0.5 - X2.z @@ -301,16 +303,16 @@ def protonateSP2(list): H = make_new_H(X2,x,y,z) H.name = "H" - + return H def make_new_H(atom, x,y,z): - new_H = Source.atom.Atom() - new_H.setProperty(numb = None, - name = 'H%s'%atom.name[1:], - resName = atom.resName, + new_H = propka.atom.Atom() + new_H.setProperty(numb = None, + name = 'H%s'%atom.name[1:], + resName = atom.resName, chainID = atom.chainID, resNumb = atom.resNumb, x = x, @@ -329,7 +331,7 @@ def make_new_H(atom, x,y,z): new_H.number_of_pi_electrons_in_double_and_triple_bonds = 0 atom.bonded_atoms.append(new_H) - atom.conformation_container.add_atom(new_H) + atom.conformation_container.add_atom(new_H) return new_H @@ -356,7 +358,7 @@ def radial_volume_desolvation(parameters, group): # ignore atoms in the same residue if atom.resNumb == group.atom.resNumb and atom.chainID == group.atom.chainID: continue - + sq_dist = squared_distance(group, atom) # desolvation @@ -375,7 +377,7 @@ def radial_volume_desolvation(parameters, group): # buried if sq_dist < parameters.buried_cutoff_squared: group.Nmass += 1 - + group.buried = calculate_weight(parameters, group.Nmass) scale_factor = calculate_scale_factor(parameters, group.buried) volume_after_allowance = max(0.00, volume-parameters.desolvationAllowance) @@ -385,9 +387,9 @@ def radial_volume_desolvation(parameters, group): # Elocl, Nlocl -> reorganisation energy (count backbone hydorgen bond acceptors, C=O) - + #print('%s %5.2f %5.2f %4d'%(group, group.buried, group.Emass, group.Nmass)) - return + return @@ -395,7 +397,7 @@ def contactDesolvation(parameters, group): """ calculates the desolvation according to the Contact Model, the old default """ - + local_radius = {'ASP': 4.5, 'GLU': 4.5, 'HIS': 4.5, @@ -405,7 +407,7 @@ def contactDesolvation(parameters, group): 'ARG': 5.0, 'C-': 4.5, '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] @@ -434,7 +436,7 @@ def contactDesolvation(parameters, group): # Note, there will be an unforseen problem: e.g. if one residue has Nmass > Nmax and # the other Nmass < Nmax, the Npair will not be Nmass1 + Nmass2! residue.buried = calculateWeight(residue.Nmass) - + return 0.00, 0.00, 0.00, 0.00 @@ -455,10 +457,10 @@ def calculate_weight(parameters, Nmass): weight = float(Nmass - parameters.Nmin)/float(parameters.Nmax - parameters.Nmin) weight = min(1.0, weight) weight = max(0.0, weight) - + return weight - + def squared_distance(atom1, atom2): # if atom1 in atom2.squared_distances: @@ -518,9 +520,9 @@ def HydrogenBondEnergy(distance, dpka_max, cutoff, f_angle=1.0): value = 0.00 else: value = 1.0-(distance-cutoff[0])/(cutoff[1]-cutoff[0]) - + dpKa = dpka_max*value*f_angle - + return abs(dpKa) @@ -535,13 +537,13 @@ def AngleFactorX(atom1=None, atom2=None, atom3=None, center=None): dX_32 = atom2.x - atom3.x dY_32 = atom2.y - atom3.y dZ_32 = atom2.z - atom3.z - + distance_23 = math.sqrt( dX_32*dX_32 + dY_32*dY_32 + dZ_32*dZ_32 ) - + dX_32 = dX_32/distance_23 dY_32 = dY_32/distance_23 dZ_32 = dZ_32/distance_23 - + if atom1 == None: dX_21 = center[0] - atom2.x dY_21 = center[1] - atom2.y @@ -550,13 +552,13 @@ def AngleFactorX(atom1=None, atom2=None, atom3=None, center=None): dX_21 = atom1.x - atom2.x dY_21 = atom1.y - atom2.y dZ_21 = atom1.z - atom2.z - + distance_12 = math.sqrt( dX_21*dX_21 + dY_21*dY_21 + dZ_21*dZ_21 ) - + dX_21 = dX_21/distance_12 dY_21 = dY_21/distance_12 dZ_21 = dZ_21/distance_12 - + f_angle = dX_21*dX_32 + dY_21*dY_32 + dZ_21*dZ_32 @@ -569,7 +571,7 @@ def hydrogen_bond_interaction(group1, group2, version): # find the smallest distance between interacting atoms atoms1 = group1.get_interaction_atoms(group2) atoms2 = group2.get_interaction_atoms(group1) - [closest_atom1, distance, closest_atom2] = Source.calculations.get_smallest_distance(atoms1, atoms2) + [closest_atom1, distance, closest_atom2] = propka.calculations.get_smallest_distance(atoms1, atoms2) if None in [closest_atom1, closest_atom2]: print('Warning: Side chain interaction failed for %s and %s'%(group1.label, group2.label)) @@ -577,20 +579,20 @@ def hydrogen_bond_interaction(group1, group2, version): # get the parameters [dpka_max, cutoff] = version.get_hydrogen_bond_parameters(closest_atom1,closest_atom2) - + if dpka_max==None or None in cutoff: return None - + # check that the closest atoms are close enough if distance >= cutoff[1]: return None # check that bond distance criteria is met - bond_distance_too_short = group1.atom.is_atom_within_bond_distance(group2.atom, + bond_distance_too_short = group1.atom.is_atom_within_bond_distance(group2.atom, version.parameters.min_bond_distance_for_hydrogen_bonds,1) if bond_distance_too_short: return None - + # set the angle factor # # ---closest_atom1/2 @@ -602,18 +604,18 @@ def hydrogen_bond_interaction(group1, group2, version): if closest_atom2.element == 'H': heavy_atom = closest_atom2.bonded_atoms[0] hydrogen = closest_atom2 - distance, f_angle, nada = Source.calculations.AngleFactorX(closest_atom1, hydrogen, heavy_atom) + distance, f_angle, nada = propka.calculations.AngleFactorX(closest_atom1, hydrogen, heavy_atom) else: # Either the structure is corrupt (no hydrogen), or the heavy atom is closer to # the titratable atom than the hydrogen. In either case we set the angle factor # to 0 f_angle = 0.0 - + elif group1.type in version.parameters.angular_dependent_sidechain_interactions: if closest_atom1.element == 'H': heavy_atom = closest_atom1.bonded_atoms[0] hydrogen = closest_atom1 - distance, f_angle, nada = Source.calculations.AngleFactorX(closest_atom2, hydrogen, heavy_atom) + distance, f_angle, nada = propka.calculations.AngleFactorX(closest_atom2, hydrogen, heavy_atom) else: # Either the structure is corrupt (no hydrogen), or the heavy atom is closer to # the titratable atom than the hydrogen. In either case we set the angle factor @@ -621,7 +623,7 @@ def hydrogen_bond_interaction(group1, group2, version): f_angle = 0.0 weight = version.calculatePairWeight(group1.Nmass, group2.Nmass) - + exception, value = version.checkExceptions(group1, group2) #exception = False # circumventing exception if exception == True: @@ -629,7 +631,7 @@ def hydrogen_bond_interaction(group1, group2, version): #print(" exception for %s %s %6.2lf" % (group1.label, group2.label, value)) else: value = version.calculateSideChainEnergy(distance, dpka_max, cutoff, weight, f_angle) - + # print('distance',distance) # print('dpka_max',dpka_max) # print('cutoff',cutoff) @@ -654,7 +656,7 @@ def HydrogenBondEnergy(distance, dpka_max, cutoff, f_angle=1.0): value = 1.0-(distance-cutoff[0])/(cutoff[1]-cutoff[0]) dpKa = dpka_max*value*f_angle - + return abs(dpKa) @@ -678,15 +680,15 @@ def checkCoulombPair(parameters, group1, group2, distance): """ Npair = group1.Nmass + group2.Nmass do_coulomb = True - + # check if both groups are titratable (ions are taken care of in determinants::setIonDeterminants) if not (group1.titratable and group2.titratable): do_coulomb = False - + # check if the distance is not too big if distance > parameters.coulomb_cutoff2: do_coulomb = False - + # check that desolvation is ok if Npair < parameters.Nmin: do_coulomb = False @@ -725,8 +727,8 @@ def BackBoneReorganization(parameters, conformation): dpKa = 0.00 for BBC_group in BBC_groups: center = [titratable_group.x, titratable_group.y, titratable_group.z] - distance, f_angle, nada = AngleFactorX(atom2=BBC_group.get_interaction_atoms(titratable_group)[0], - atom3=BBC_group.atom, + distance, f_angle, nada = AngleFactorX(atom2=BBC_group.get_interaction_atoms(titratable_group)[0], + atom3=BBC_group.atom, center=center) if distance < 6.0 and f_angle > 0.001: value = 1.0-(distance-3.0)/(6.0-3.0) @@ -767,7 +769,7 @@ def checkExceptions(version, group1, group2): else: # do nothing, no exception for this pair exception = False; value = None - + return exception, value @@ -780,7 +782,7 @@ def checkCooArgException(group_coo, group_arg, version): str = "xxx" exception = True value_tot = 0.00 - + #dpka_max = parameters.sidechain_interaction.get_value(group_coo.type,group_arg.type) #cutoff = parameters.sidechain_cutoffs.get_value(group_coo.type,group_arg.type) @@ -801,7 +803,7 @@ def checkCooArgException(group_coo, group_arg, version): if group_arg.type in version.parameters.angular_dependent_sidechain_interactions: atom3 = closest_arg_atom.bonded_atoms[0] distance, f_angle, nada = AngleFactorX(closest_coo_atom, closest_arg_atom, atom3) - + value = HydrogenBondEnergy(distance, dpka_max, cutoff, f_angle) #print(iter, closest_coo_atom, closest_arg_atom,distance,value) value_tot += value @@ -809,7 +811,7 @@ def checkCooArgException(group_coo, group_arg, version): atoms_coo.remove(closest_coo_atom) atoms_arg.remove(closest_arg_atom) - + return exception, value_tot @@ -820,7 +822,7 @@ def checkCooCooException(group1, group2, version): exception = True [closest_atom1, distance, closest_atom2] = get_smallest_distance(group1.get_interaction_atoms(group2), group2.get_interaction_atoms(group1)) - + #dpka_max = parameters.sidechain_interaction.get_value(group1.type,group2.type) #cutoff = parameters.sidechain_cutoffs.get_value(group1.type,group2.type) [dpka_max, cutoff] = version.get_hydrogen_bond_parameters(closest_atom1,closest_atom2) @@ -863,7 +865,7 @@ def checkCysHisException(group1, group2, version): if checkBuried(group1.Nmass, group2.Nmass): exception = True - return exception, version.parameters.CYS_HIS_exception + return exception, version.parameters.CYS_HIS_exception def checkCysCysException(group1, group2, version): diff --git a/Source/conformation_container.py b/propka/conformation_container.py similarity index 92% rename from Source/conformation_container.py rename to propka/conformation_container.py index c327390..76439d4 100644 --- a/Source/conformation_container.py +++ b/propka/conformation_container.py @@ -2,7 +2,10 @@ # Container for molecular conformations # -import Source.group, Source.determinants, Source.determinant, Source.ligand, Source.output, Source.coupled_groups, functools +from __future__ import division +from __future__ import print_function + +import propka.group, propka.determinants, propka.determinant, propka.ligand, propka.output, propka.coupled_groups, functools class Conformation_container: def __init__(self, name='', parameters=None, molecular_container=None): @@ -13,10 +16,10 @@ class Conformation_container: self.groups = [] self.chains = [] self.current_iter_item = 0 - + self.marvin_pkas_calculated = False self.non_covalently_coupled_groups = False - + return @@ -28,21 +31,21 @@ class Conformation_container: for atom in self.get_non_hydrogen_atoms(): # has this atom been checked for groups? if atom.groups_extracted == 0: - self.validate_group(Source.group.is_group(self.parameters, atom)) - # if the atom has been checked in a another conformation, check if it has a - # group that should be used in this conformation as well + self.validate_group(propka.group.is_group(self.parameters, atom)) + # if the atom has been checked in a another conformation, check if it has a + # group that should be used in this conformation as well elif atom.group: self.validate_group(atom.group) - return + return def additional_setup_when_reading_input_file(self): - - # if a group is coupled and we are reading a .propka_input file, + + # if a group is coupled and we are reading a .propka_input file, # some more configuration might be needed # print coupling map - map = Source.output.make_interaction_map('Covalent coupling map for %s'%self, + map = propka.output.make_interaction_map('Covalent coupling map for %s'%self, self.get_covalently_coupled_groups(), lambda g1,g2: g1 in g2.covalently_coupled_groups) print(map) @@ -50,11 +53,11 @@ class Conformation_container: # check if we should set a common charge centre as well if self.parameters.common_charge_centre: self.set_common_charge_centres() - + return - + def set_common_charge_centres(self): - for system in self.get_coupled_systems(self.get_covalently_coupled_groups(), Source.group.Group.get_covalently_coupled_groups): + for system in self.get_coupled_systems(self.get_covalently_coupled_groups(), propka.group.Group.get_covalently_coupled_groups): # make a list of the charge centre coordinates all_coordinates = list(map(lambda g: [g.x, g.y, g.z], system)) # find the common charge center @@ -82,19 +85,19 @@ class Conformation_container: continue if cg.atom.sybyl_type == group.atom.sybyl_type: group.couple_covalently(cg) - + # check if we should set a common charge centre as well if self.parameters.common_charge_centre: self.set_common_charge_centres() # print coupling map - map = Source.output.make_interaction_map('Covalent coupling map for %s'%self, + map = propka.output.make_interaction_map('Covalent coupling map for %s'%self, #self.get_titratable_groups(), self.get_covalently_coupled_groups(), lambda g1,g2: g1 in g2.covalently_coupled_groups) print(map) - + return @@ -102,8 +105,8 @@ class Conformation_container: # check if non-covalent coupling has already been set up in an input file if len(list(filter(lambda g: len(g.non_covalently_coupled_groups)>0, self.get_titratable_groups())))>0: self.non_covalently_coupled_groups = True - - Source.coupled_groups.nccg.identify_non_covalently_coupled_groups(self,verbose=verbose) + + propka.coupled_groups.nccg.identify_non_covalently_coupled_groups(self,verbose=verbose) # re-do the check if len(list(filter(lambda g: len(g.non_covalently_coupled_groups)>0, self.get_titratable_groups())))>0: @@ -157,26 +160,26 @@ class Conformation_container: def calculate_pka(self, version, options): print('\nCalculating pKas for',self) - + # calculate desolvation for group in self.get_titratable_groups()+self.get_ions(): version.calculate_desolvation(group) # calculate backbone interactions - Source.determinants.setBackBoneDeterminants(self.get_titratable_groups(), self.get_backbone_groups(), version) + propka.determinants.setBackBoneDeterminants(self.get_titratable_groups(), self.get_backbone_groups(), version) # setting ion determinants - Source.determinants.setIonDeterminants(self, version) + propka.determinants.setIonDeterminants(self, version) # calculating the back-bone reorganization/desolvation term version.calculateBackBoneReorganization(self) # setting remaining non-iterative and iterative side-chain & Coulomb interaction determinants - Source.determinants.setDeterminants(self.get_sidechain_groups(), version=version, options=options) + propka.determinants.setDeterminants(self.get_sidechain_groups(), version=version, options=options) # calculating the total pKa values for group in self.groups: group.calculate_total_pka() - + # take coupling effects into account penalised_labels = self.coupling_effects() @@ -195,20 +198,20 @@ class Conformation_container: def coupling_effects(self): # - # Bases: The group with the highest pKa (the most stable one in the - # charged form) will be the first one to adopt a proton as pH + # Bases: The group with the highest pKa (the most stable one in the + # charged form) will be the first one to adopt a proton as pH # is lowered and this group is allowed to titrate. # The remaining groups are penalised # # Acids: The group with the highest pKa (the least stable one in the - # charged form) will be the last group to loose the proton as + # charged form) will be the last group to loose the proton as # pH is raised and will be penalised. # The remaining groups are allowed to titrate. # penalised_labels = [] - - for all_groups in self.get_coupled_systems(self.get_covalently_coupled_groups(), - Source.group.Group.get_covalently_coupled_groups): + + for all_groups in self.get_coupled_systems(self.get_covalently_coupled_groups(), + propka.group.Group.get_covalently_coupled_groups): # check if we should share determinants if self.parameters.shared_determinants: @@ -227,14 +230,14 @@ class Conformation_container: for a in all_groups: if a == first_group: continue # group with the highest pKa is allowed to titrate... - a.coupled_titrating_group = first_group + a.coupled_titrating_group = first_group penalised_labels.append(a.label) #... and the rest is penalised return penalised_labels def share_determinants(self, groups): - + # make a list of the determinants to share types = ['sidechain','backbone','coulomb'] for type in types: @@ -250,7 +253,7 @@ class Conformation_container: # overwrite/add maximum value for each determinant for det_group in max_dets.keys(): - new_determinant = Source.determinant.Determinant(det_group, max_dets[det_group]) + new_determinant = propka.determinant.Determinant(det_group, max_dets[det_group]) for g in groups: g.set_determinant(new_determinant,type) @@ -260,7 +263,7 @@ class Conformation_container: def get_coupled_systems(self, groups, get_coupled_groups): """ This generator will yield one covalently coupled system at the time """ - groups = set(groups) + groups = set(groups) while len(groups)>0: # extract a system of coupled groups ... system = set() @@ -272,7 +275,7 @@ class Conformation_container: return - + def get_a_coupled_system_of_groups(self, new_group, coupled_groups, get_coupled_groups): coupled_groups.add(new_group) [self.get_a_coupled_system_of_groups(c, coupled_groups, get_coupled_groups) for c in get_coupled_groups(new_group) if c not in coupled_groups] @@ -311,7 +314,7 @@ class Conformation_container: """ returns all sidechain groups needed for the pKa calculations """ return [group for group in self.groups if ('BB' not in group.type\ and not group.atom.cysteine_bridge)] - + def get_covalently_coupled_groups(self): return [g for g in self.groups if len(g.covalently_coupled_groups)>0] @@ -331,7 +334,7 @@ class Conformation_container: def get_titratable_groups(self): return [group for group in self.groups if group.titratable] - + def get_titratable_groups_and_cysteine_bridges(self): return [group for group in self.groups if group.titratable or group.residue_type == 'CYS'] @@ -367,19 +370,19 @@ class Conformation_container: atom.conformation_container = self if not atom.molecular_container: atom.molecular_container = self.molecular_container - + # store chain id for bookkeeping if not atom.chainID in self.chains: self.chains.append(atom.chainID) - - return + + return def copy_atom(self, atom): new_atom = atom.makeCopy() self.atoms.append(new_atom) new_atom.conformation_container = self - return + return def get_non_hydrogen_atoms(self): return [atom for atom in self.atoms if atom.element!='H'] @@ -391,7 +394,7 @@ class Conformation_container: if not self.have_atom(atom): self.copy_atom(atom) return - + def have_atom(self, atom): res = False for a in self.atoms: @@ -407,12 +410,12 @@ class Conformation_container: return g return False - + def set_ligand_atom_names(self): for atom in self.get_ligand_atoms(): - Source.ligand.assign_sybyl_type(atom) + propka.ligand.assign_sybyl_type(atom) return - + def __str__(self): @@ -431,7 +434,7 @@ class Conformation_container: return def sort_atoms_key(self, atom): - key = ord(atom.chainID)*1e7 + key = ord(atom.chainID)*1e7 key += atom.resNumb*1000 if len(atom.name) > len(atom.element): key += ord(atom.name[len(atom.element)]) diff --git a/Source/coupled_groups.py b/propka/coupled_groups.py similarity index 95% rename from Source/coupled_groups.py rename to propka/coupled_groups.py index 06d966f..1c2fa7b 100644 --- a/Source/coupled_groups.py +++ b/propka/coupled_groups.py @@ -1,12 +1,15 @@ -import math, Source.output, Source.group, Source.lib, itertools +from __future__ import division +from __future__ import print_function + +import math, propka.output, propka.group, propka.lib, itertools class non_covalently_couple_groups: def __init__(self): self.parameters = None - + # self.do_intrinsic = False # self.do_pair_wise = False self.do_prot_stat = True @@ -18,10 +21,10 @@ class non_covalently_couple_groups: # Methods for finding coupled groups # def is_coupled_protonation_state_probability(self, group1, group2, energy_method, return_on_fail=True): - + # check if the interaction energy is high enough interaction_energy = max(self.get_interaction(group1,group2), self.get_interaction(group2,group1)) - + if interaction_energy<=self.parameters.min_interaction_energy and return_on_fail: return {'coupling_factor':-1.0} @@ -29,7 +32,7 @@ class non_covalently_couple_groups: for group in [group1, group2]: if not hasattr(group, 'intrinsic_pKa'): group.calculate_intrinsic_pka() - + use_pH = self.parameters.pH if self.parameters.pH == 'variable': use_pH = min(group1.pka_value, group2.pka_value) @@ -61,7 +64,7 @@ class non_covalently_couple_groups: self.swap_interactions([group1], [group2]) group1.calculate_total_pka() group2.calculate_total_pka() - + # check difference in free energy if abs(default_energy - swapped_energy) > self.parameters.max_free_energy_diff and return_on_fail: return {'coupling_factor':-1.0} @@ -80,21 +83,21 @@ class non_covalently_couple_groups: self.get_interaction_factor(interaction_energy) return {'coupling_factor':factor, - 'default_energy':default_energy, - 'swapped_energy':swapped_energy, + 'default_energy':default_energy, + 'swapped_energy':swapped_energy, 'interaction_energy':interaction_energy, - 'swapped_pka1':swapped_pka1, + 'swapped_pka1':swapped_pka1, 'swapped_pka2':swapped_pka2, 'pka_shift1':pka_shift1, 'pka_shift2':pka_shift2, 'pH':use_pH} - + # # Methods for calculating the coupling factor - # - + # + def get_pka_diff_factor(self, pka1, pka2): intrinsic_pka_diff = abs(pka1-pka2) res = 0.0 @@ -120,10 +123,10 @@ class non_covalently_couple_groups: - + def identify_non_covalently_coupled_groups(self, conformation, verbose=True): """ Finds coupled residues in protein """ - + self.parameters = conformation.parameters if verbose: print('') @@ -156,19 +159,19 @@ class non_covalently_couple_groups: data = self.is_coupled_protonation_state_probability(g1, g2,conformation.calculate_folding_energy) if data['coupling_factor'] >0.0: g1.couple_non_covalently(g2) - + if verbose: self.print_out_swaps(conformation) return def print_out_swaps(self, conformation, verbose=True): - map = Source.output.make_interaction_map('Non-covalent coupling map for %s'%conformation, + map = propka.output.make_interaction_map('Non-covalent coupling map for %s'%conformation, conformation.get_non_covalently_coupled_groups(), lambda g1,g2: g1 in g2.non_covalently_coupled_groups) print(map) - for system in conformation.get_coupled_systems(conformation.get_non_covalently_coupled_groups(),Source.group.Group.get_non_covalently_coupled_groups): + for system in conformation.get_coupled_systems(conformation.get_non_covalently_coupled_groups(),propka.group.Group.get_non_covalently_coupled_groups): self.print_system(conformation, list(system)) return @@ -187,7 +190,7 @@ class non_covalently_couple_groups: print(coup_info) # make list of possible combinations of swap to try out - combinations = Source.lib.generate_combinations(interactions) + combinations = propka.lib.generate_combinations(interactions) # Make possible swap combinations swap_info = '' @@ -199,7 +202,7 @@ class non_covalently_couple_groups: for interaction in combination: swap_info += ' %s %s\n'%(interaction[0].label, interaction[1].label) - # swap... + # swap... for interaction in combination: self.swap_interactions([interaction[0]],[interaction[1]]) @@ -213,13 +216,13 @@ class non_covalently_couple_groups: return # # Interaction and swapping methods - # + # def get_interaction(self, group1, group2, include_side_chain_hbs = True): determinants = group1.determinants['coulomb'] if include_side_chain_hbs: determinants = group1.determinants['sidechain'] + group1.determinants['coulomb'] - + interaction_energy = 0.0 for det in determinants: if group2 == det.group: @@ -234,7 +237,7 @@ class non_covalently_couple_groups: s = ' '+'-'*113+'\n' for group in system: s += self.tagged_print(' %-8s|'%tag,group.getDeterminantString(), all_labels) - + return s+'\n' def swap_interactions(self, groups1, groups2, include_side_chain_hbs = True): @@ -242,7 +245,7 @@ class non_covalently_couple_groups: for i in range(len(groups1)): group1 = groups1[i] group2 = groups2[i] - + # swap the interactions! self.transfer_determinant(group1.determinants['coulomb'], group2.determinants['coulomb'], group1.label, group2.label) if include_side_chain_hbs: @@ -266,7 +269,7 @@ class non_covalently_couple_groups: for det in determinants2: if det.label == label1: from2to1.append(det) - + # ...and transfer it! for det in from1to2: det.label = label1 @@ -280,7 +283,7 @@ class non_covalently_couple_groups: return - # + # # Output methods # @@ -295,7 +298,7 @@ class non_covalently_couple_groups: def make_data_to_string(self, data, group1, group2): s = """ %s and %s coupled (prot.state): %5.2f Energy levels: %6.2f, %6.2f (difference: %6.2f) at pH %6.2f - Interaction energy: %6.2f + Interaction energy: %6.2f Intrinsic pka's: %6.2f, %6.2f (difference: %6.2f) Swapped pKa's: %6.2f, %6.2f (difference: %6.2f, %6.2f)"""%(group1.label, group2.label, diff --git a/Source/determinant.py b/propka/determinant.py similarity index 88% rename from Source/determinant.py rename to propka/determinant.py index 3be7a02..5085ca1 100644 --- a/Source/determinant.py +++ b/propka/determinant.py @@ -1,4 +1,7 @@ +from __future__ import division +from __future__ import print_function + class Determinant: """ Determinant class - set up for later structurization diff --git a/Source/determinants.py b/propka/determinants.py similarity index 91% rename from Source/determinants.py rename to propka/determinants.py index 8f2dc60..b3ed586 100644 --- a/Source/determinants.py +++ b/propka/determinants.py @@ -1,8 +1,12 @@ + +from __future__ import division +from __future__ import print_function + import math, time -import Source.iterative, Source.lib, Source.vector_algebra -import Source.calculations -from Source.determinant import Determinant +import propka.iterative, propka.lib, propka.vector_algebra +import propka.calculations +from propka.determinant import Determinant def setDeterminants(propka_groups, version=None, options=None): @@ -19,19 +23,19 @@ def setDeterminants(propka_groups, version=None, options=None): # do not calculate interactions for coupled groups if group2 in group1.covalently_coupled_groups: break - - distance = Source.calculations.distance(group1, group2) - + + distance = propka.calculations.distance(group1, group2) + if distance < version.parameters.coulomb_cutoff2: interaction_type = version.parameters.interaction_matrix.get_value(group1.type,group2.type) if interaction_type == 'I': - Source.iterative.addtoDeterminantList(group1, group2, distance, iterative_interactions, version=version) + propka.iterative.addtoDeterminantList(group1, group2, distance, iterative_interactions, version=version) elif interaction_type == 'N': addDeterminants(group1, group2, distance, version) - + # --- Iterative section ---# - Source.iterative.addDeterminants(iterative_interactions, version, options=options) + propka.iterative.addDeterminants(iterative_interactions, version, options=options) def addDeterminants(group1, group2, distance, version): @@ -52,11 +56,11 @@ def addSidechainDeterminants(group1, group2, version=None): adding side-chain determinants/perturbations Note, resNumb1 > resNumb2 """ - + hbond_interaction = version.hydrogen_bond_interaction(group1, group2) if hbond_interaction: - + if group1.charge == group2.charge: # acid pair or base pair if group1.model_pka < group2.model_pka: @@ -78,7 +82,7 @@ def addCoulombDeterminants(group1, group2, distance, version): """ adding NonIterative Coulomb determinants/perturbations """ - + coulomb_interaction = version.electrostatic_interaction(group1, group2, distance) if coulomb_interaction: @@ -104,7 +108,7 @@ def addCoulombAcidPair(object1, object2, value): Adding the Coulomb interaction (an acid pair): the higher pKa is raised """ - + if object1.model_pka > object2.model_pka: newDeterminant = Determinant(object2, value) object1.determinants['coulomb'].append(newDeterminant) @@ -151,7 +155,7 @@ def setIonDeterminants(conformation_container, version): """ for titratable_group in conformation_container.get_titratable_groups(): for ion_group in conformation_container.get_ions(): - squared_distance = Source.calculations.squared_distance(titratable_group, ion_group) + squared_distance = propka.calculations.squared_distance(titratable_group, ion_group) if squared_distance < version.parameters.coulomb_cutoff2_squared: weight = version.calculatePairWeight(titratable_group.Nmass, ion_group.Nmass) # the pKa of both acids and bases are shifted up by negative ions (and vice versa) @@ -162,16 +166,16 @@ def setIonDeterminants(conformation_container, version): return def setBackBoneDeterminants(titratable_groups, backbone_groups, version): - + for titratable_group in titratable_groups: # find out which backbone groups this titratable is interacting with for backbone_group in backbone_groups: # find the interacting atoms backbone_interaction_atoms = backbone_group.get_interaction_atoms(titratable_group) titratable_group_interaction_atoms = titratable_group.interaction_atoms_for_acids - + # find the smallest distance - [backbone_atom, distance, titratable_atom] = Source.calculations.get_smallest_distance(backbone_interaction_atoms, + [backbone_atom, distance, titratable_atom] = propka.calculations.get_smallest_distance(backbone_interaction_atoms, titratable_group_interaction_atoms) # get the parameters parameters = version.get_backbone_hydrogen_bond_parameters(backbone_atom, titratable_atom) @@ -186,7 +190,7 @@ def setBackBoneDeterminants(titratable_groups, backbone_groups, version): # for BBC groups, the hydrogen is on the titratable group # # Titra. - # / + # / # H # . # O @@ -197,7 +201,7 @@ def setBackBoneDeterminants(titratable_groups, backbone_groups, version): if titratable_atom.element == 'H': heavy_atom = titratable_atom.bonded_atoms[0] hydrogen_atom = titratable_atom - [d1, f_angle, d2] = Source.calculations.AngleFactorX(atom1=heavy_atom, + [d1, f_angle, d2] = propka.calculations.AngleFactorX(atom1=heavy_atom, atom2=hydrogen_atom, atom3=backbone_atom) else: @@ -218,7 +222,7 @@ def setBackBoneDeterminants(titratable_groups, backbone_groups, version): if backbone_atom.element == 'H': backbone_N = backbone_atom.bonded_atoms[0] backbone_H = backbone_atom - [d1, f_angle, d2] = Source.calculations.AngleFactorX(atom1=titratable_atom, + [d1, f_angle, d2] = propka.calculations.AngleFactorX(atom1=titratable_atom, atom2=backbone_H, atom3=backbone_N) else: @@ -226,13 +230,13 @@ def setBackBoneDeterminants(titratable_groups, backbone_groups, version): # the titratable atom than the hydrogen. In either case we set the angle factor # to 0 f_angle = 0.0 - - + + if f_angle > 0.001: - value = titratable_group.charge * Source.calculations.HydrogenBondEnergy(distance, dpKa_max, [cutoff1,cutoff2], f_angle) - + value = titratable_group.charge * propka.calculations.HydrogenBondEnergy(distance, dpKa_max, [cutoff1,cutoff2], f_angle) + newDeterminant = Determinant(backbone_group, value) titratable_group.determinants['backbone'].append(newDeterminant) - + return diff --git a/Source/group.py b/propka/group.py similarity index 96% rename from Source/group.py rename to propka/group.py index 105b742..ac52d3e 100644 --- a/Source/group.py +++ b/propka/group.py @@ -2,9 +2,12 @@ # Class for storing groups important for propka calculations # -import Source.ligand, Source.determinant, Source.ligand_pka_values, math, Source.protonate +from __future__ import division +from __future__ import print_function -my_protonator = Source.protonate.Protonate(verbose=False) +import propka.ligand, propka.determinant, propka.ligand_pka_values, math, propka.protonate + +my_protonator = propka.protonate.Protonate(verbose=False) expected_atoms_acid_interactions = { 'COO':{'O':2}, @@ -65,7 +68,7 @@ expected_atoms_base_interactions = { 'SH': {'S':1}, 'CG': {'N':3}, 'C2N':{'N':2}, - 'OCO':{'O':2}, + 'OCO':{'O':2}, 'N30':{'N':1}, 'N31':{'N':1}, 'N32':{'N':1}, @@ -86,7 +89,7 @@ class Group: self.determinants = {'sidechain':[],'backbone':[],'coulomb':[]} self.pka_value = 0.0 self.model_pka = 0.0 - + self.Emass = 0.0 self.Nmass = 0.0 self.Elocl = 0.0 @@ -99,18 +102,18 @@ class Group: self.interaction_atoms_for_acids = [] self.interaction_atoms_for_bases = [] self.model_pka_set = False - + # information on covalent and non-covalent coupling self.non_covalently_coupled_groups = [] self.covalently_coupled_groups = [] self.coupled_titrating_group = None self.common_charge_centre = False - + self.residue_type = self.atom.resName 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 ']: @@ -119,7 +122,7 @@ class Group: atom.name.replace('\'','')[-1], atom.resNumb, atom.chainID) - + # self.label = '%1s%1s%1s%4d%2s'%(self.residue_type[1], atom.element,atom.name[-1], atom.resNumb, atom.chainID) else: self.label = '%-3s%4s%2s'%(self.residue_type, atom.name, atom.chainID) @@ -133,28 +136,28 @@ class Group: # # Coupling-related methods - # + # def couple_covalently(self, other): """ Couple this group with another group """ # do the coupling if not other in self.covalently_coupled_groups: self.covalently_coupled_groups.append(other) - + if not self in other.covalently_coupled_groups: other.covalently_coupled_groups.append(self) - return + return def couple_non_covalently(self, other): """ Couple this group with another group """ # do the coupling if not other in self.non_covalently_coupled_groups: self.non_covalently_coupled_groups.append(other) - + if not self in other.non_covalently_coupled_groups: other.non_covalently_coupled_groups.append(self) - return + return def get_covalently_coupled_groups(self): return self.covalently_coupled_groups def get_non_covalently_coupled_groups(self): return self.non_covalently_coupled_groups @@ -168,12 +171,12 @@ class Group: continue for type in ['sidechain','backbone','coulomb']: - for g in other.determinants[type]: self.share_determinant(g,type) - + for g in other.determinants[type]: self.share_determinant(g,type) + # recalculate pka values self.calculate_total_pka() other.calculate_total_pka() - + return @@ -190,19 +193,19 @@ class Group: # otherwise we just add the determinant to our list if not added: - self.determinants[type].append(Source.determinant.Determinant(new_determinant.group, + self.determinants[type].append(propka.determinant.Determinant(new_determinant.group, new_determinant.value)) return def make_covalently_coupled_line(self): - + # first check if there are any coupled groups at all if len(self.covalently_coupled_groups) == 0: return '' line = 'CCOUPL%5d'%self.atom.numb - + # extract and sort numbers of coupled groups coupled = [] for g in self.covalently_coupled_groups: @@ -213,7 +216,7 @@ class Group: for b in coupled: line += '%5d'%b line += '\n' - + return line def make_non_covalently_coupled_line(self): @@ -222,7 +225,7 @@ class Group: return '' line = 'NCOUPL%5d'%self.atom.numb - + # extract and sort numbers of coupled groups coupled = [] for g in self.non_covalently_coupled_groups: @@ -233,7 +236,7 @@ class Group: for b in coupled: line += '%5d'%b line += '\n' - + return line @@ -283,7 +286,7 @@ class Group: return # otherwise we just add the determinant to our list - self.determinants[type].append(Source.determinant.Determinant(new_determinant.group, + self.determinants[type].append(propka.determinant.Determinant(new_determinant.group, new_determinant.value)) return @@ -298,7 +301,7 @@ class Group: return # otherwise we just add the determinant to our list - self.determinants[type].append(Source.determinant.Determinant(new_determinant.group, + self.determinants[type].append(propka.determinant.Determinant(new_determinant.group, new_determinant.value)) return @@ -344,12 +347,12 @@ class Group: self.charge = self.parameters.charge[self.type] if self.residue_type in self.parameters.ions.keys(): self.charge = self.parameters.ions[self.residue_type] - + #print('ION setup',self,self.residue_type, self.charge) # find the center and the interaction atoms self.setup_atoms() - + # set the model pka value self.titratable = False if self.residue_type in self.parameters.model_pkas.keys(): @@ -365,7 +368,7 @@ class Group: if self.model_pka_set and not self.atom.cysteine_bridge: self.titratable = True - + return def setup_atoms(self): @@ -378,7 +381,7 @@ class Group: def set_interaction_atoms(self, interaction_atoms_for_acids, interaction_atoms_for_bases): [a.set_group_type(self.type) for a in interaction_atoms_for_acids+interaction_atoms_for_bases] - + self.interaction_atoms_for_acids = interaction_atoms_for_acids self.interaction_atoms_for_bases = interaction_atoms_for_bases @@ -400,11 +403,11 @@ class Group: print(' Expected %d interaction atoms for acids, found:'%Na) for i in range(len(self.interaction_atoms_for_acids)): print(' %s'%self.interaction_atoms_for_acids[i]) - + print(' Expected %d interaction atoms for bases, found:'%Nb) for i in range(len(self.interaction_atoms_for_bases)): print(' %s'%self.interaction_atoms_for_bases[i]) - + #return @@ -437,7 +440,7 @@ class Group: number_of_sidechain = len(self.determinants['sidechain']) number_of_backbone = len(self.determinants['backbone']) number_of_coulomb = len(self.determinants['coulomb']) - + number_of_lines = max(1, number_of_sidechain, number_of_backbone, number_of_coulomb) str = "" for line_number in range(number_of_lines): @@ -458,11 +461,11 @@ class Group: str += " %6.2lf %4d" % (self.Elocl, self.Nlocl) else: str += "%40s" % (" ") - + # add the determinants for type in ['sidechain','backbone','coulomb']: str += self.get_determinant_for_string(type,line_number) - + # adding end-of-line str += "\n" @@ -484,7 +487,7 @@ class Group: if self.atom.cysteine_bridge: self.pka_value = 99.99 return - + self.pka_value = self.model_pka + self.Emass + self.Elocl @@ -494,12 +497,12 @@ class Group: return - + def calculate_intrinsic_pka(self): """ - Calculates the intrinsic pKa values from the desolvation determinants, back-bone hydrogen bonds, + 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.0 @@ -528,18 +531,18 @@ class Group: ligand_type = '' if self.atom.type == 'hetatm': ligand_type = self.type - + penalty = '' if self.coupled_titrating_group: penalty = ' NB: Discarded due to coupling with %s'%self.coupled_titrating_group.label - str = " %9s %8.2lf %10.2lf %18s %s\n" % (self.label, - self.pka_value, - self.model_pka,ligand_type, + str = " %9s %8.2lf %10.2lf %18s %s\n" % (self.label, + self.pka_value, + self.model_pka,ligand_type, penalty) return str - + def __str__(self): return 'Group (%s) for %s'%(self.type,self.atom) @@ -559,7 +562,7 @@ class Group: reference = parameters.reference # If not titratable, the contribution is zero - + if not self.titratable: return 0.00 @@ -571,25 +574,25 @@ class Group: if determinant.value > 0.00: pka_prime -= determinant.value ddG_neutral = -1.36*(pka_prime - self.model_pka) - + # calculating the ddG(low-pH --> pH) contribution # folded x = pH - self.pka_value y = 10**x Q_pro = math.log10(1+y) - + # unfolded x = pH - self.model_pka 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 calculate_charge(self, parmaeters, pH=7.0, state='folded'): - + if state == "unfolded": x = self.charge * (self.model_pka - pH) else: @@ -597,7 +600,7 @@ class Group: y = 10**x charge = self.charge*(y/(1.0+y)) - + return charge @@ -610,7 +613,7 @@ class COO_group(Group): def setup_atoms(self): # Identify the two caroxyl oxygen atoms the_oxygens = self.atom.get_bonded_elements('O') - + # set the center using the two oxygen carboxyl atoms self.set_center(the_oxygens) self.set_interaction_atoms(the_oxygens, the_oxygens) @@ -624,7 +627,7 @@ class HIS_group(Group): def setup_atoms(self): # Find the atoms in the histidine ring - ring_atoms = Source.ligand.is_ring_member(self.atom) + ring_atoms = propka.ligand.is_ring_member(self.atom) if len(ring_atoms) != 5: print('Warning: His group does not seem to contain a ring',self) @@ -634,11 +637,11 @@ class HIS_group(Group): # set the center using the ring atoms self.set_center(ring_atoms) - + # find the hydrogens on the ring-nitrogens hydrogens = [] nitrogens = [ra for ra in ring_atoms if ra.element == 'N'] - + for nitrogen in nitrogens: hydrogens.extend(nitrogen.get_bonded_elements('H')) @@ -673,12 +676,12 @@ class ARG_group(Group): def setup_atoms(self): # set the center at the position of the main atom self.set_center([self.atom]) - + # find all the hydrogens on the nitrogen atoms nitrogens = self.atom.get_bonded_elements('N') for n in nitrogens: my_protonator.protonate_atom(n) - + hydrogens = [] for nitrogen in nitrogens: hydrogens.extend(nitrogen.get_bonded_elements('H')) @@ -705,19 +708,19 @@ class AMD_group(Group): def setup_atoms(self): # Identify the oxygen and nitrogen amide atoms the_oxygen = self.atom.get_bonded_elements('O') - the_nitrogen = self.atom.get_bonded_elements('N') - - # add protons to the nitrogen + the_nitrogen = self.atom.get_bonded_elements('N') + + # add protons to the nitrogen my_protonator.protonate_atom(the_nitrogen[0]) the_hydrogens = the_nitrogen[0].get_bonded_elements('H') # set the center using the oxygen and nitrogen amide atoms self.set_center(the_oxygen+the_nitrogen) - + self.set_interaction_atoms(the_nitrogen+the_hydrogens,the_oxygen) return - + class TRP_group(Group): def __init__(self, atom): @@ -727,7 +730,7 @@ class TRP_group(Group): def setup_atoms(self): # set the center at the position of the main atom self.set_center([self.atom]) - + # find the hydrogen on the nitrogen atom my_protonator.protonate_atom(self.atom) the_hydrogen = self.atom.get_bonded_elements('H') @@ -751,12 +754,12 @@ class Cterm_group(Group): the_carbon = self.atom.get_bonded_elements('C') the_other_oxygen = the_carbon[0].get_bonded_elements('O') the_other_oxygen.remove(self.atom) - + # set the center and interaction atoms the_oxygens = [self.atom]+ the_other_oxygen self.set_center(the_oxygens) self.set_interaction_atoms(the_oxygens, the_oxygens) - + return @@ -802,7 +805,7 @@ class NAR_group(Group): self.residue_type = 'NAR' print('Found NAR group:',atom) return - + def setup_atoms(self): # Identify the hydrogen @@ -825,7 +828,7 @@ class NAM_group(Group): self.residue_type = 'NAM' print('Found NAM group:',atom) return - + def setup_atoms(self): # Identify the hydrogen @@ -936,7 +939,7 @@ class CG_group(Group): # set the center using the nitrogen self.set_center([self.atom]) - + the_hydrogens = [] for n in the_nitrogens: my_protonator.protonate_atom(n) @@ -1065,7 +1068,7 @@ class NP1_group(Group): print('Found NP1 group:',atom) return - + def setup_atoms(self): # Identify the nitrogens my_protonator.protonate_atom(self.atom) @@ -1123,16 +1126,16 @@ class titratable_ligand_group(Group): self.model_pka = atom.marvin_pka print('Titratable ligand group ',atom, self.model_pka, self.charge) self.model_pka_set = True - + return def is_group(parameters, atom): atom.groups_extracted = 1 - + # check if this atom belongs to a protein group protein_group = is_protein_group(parameters, atom) - if protein_group: return protein_group + if protein_group: return protein_group # check if this atom belongs to a ion group ion_group = is_ion_group(parameters, atom) @@ -1168,7 +1171,7 @@ 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.resName != 'PRO': return BBN_group(atom) if atom.type == 'atom' and atom.name == 'C': # ignore C- carboxyl @@ -1199,10 +1202,10 @@ def is_ligand_group_by_groups(parameters, atom): if atom.sybyl_type == 'N.ar': if len(atom.get_bonded_heavy_atoms())==2: return NAR_group(atom) - + if atom.sybyl_type == 'N.am': return NAM_group(atom) - + if atom.sybyl_type in ['N.3', 'N.4']: heavy_bonded = atom.get_bonded_heavy_atoms() if len(heavy_bonded) == 0: @@ -1218,14 +1221,14 @@ def is_ligand_group_by_groups(parameters, atom): return N1_group(atom) if atom.sybyl_type == 'N.pl3': - # make sure that this atom is not part of a guadinium or amidinium group + # make sure that this atom is not part of a guadinium or amidinium group bonded_carbons = atom.get_bonded_elements('C') if len(bonded_carbons) == 1: bonded_nitrogens = bonded_carbons[0].get_bonded_elements('N') if len(bonded_nitrogens) == 1: return NP1_group(atom) - + if atom.sybyl_type == 'C.2': # Guadinium and amidinium groups bonded_nitrogens = atom.get_bonded_elements('N') @@ -1263,7 +1266,7 @@ def is_ligand_group_by_groups(parameters, atom): if atom.sybyl_type == 'O.2': return O2_group(atom) - + if atom.sybyl_type == 'S.3': # sulphide group @@ -1272,8 +1275,8 @@ def is_ligand_group_by_groups(parameters, atom): # other SP3 oxygen #else: # return S3_group(atom) - - + + return None @@ -1284,11 +1287,11 @@ def is_ligand_group_by_marvin_pkas(parameters, atom): # calculate Marvin ligand pkas for this conformation container # if not already done if not atom.conformation_container.marvin_pkas_calculated: - lpka = Source.ligand_pka_values.ligand_pka_values(parameters) + lpka = propka.ligand_pka_values.ligand_pka_values(parameters) lpka.get_marvin_pkas_for_molecular_container(atom.molecular_container, min_pH=parameters.min_ligand_model_pka, max_pH=parameters.max_ligand_model_pka) - + if atom.marvin_pka: return titratable_ligand_group(atom) diff --git a/Source/iterative.py b/propka/iterative.py similarity index 98% rename from Source/iterative.py rename to propka/iterative.py index c994867..f2316bc 100644 --- a/Source/iterative.py +++ b/propka/iterative.py @@ -1,7 +1,12 @@ + +from __future__ import division +from __future__ import print_function + import math, time -import Source.lib as lib -from Source.determinant import Determinant -import Source.calculations + +import propka.lib as lib +from propka.determinant import Determinant +import propka.calculations # Some library functions for the interative pKa determinants @@ -32,11 +37,11 @@ def addtoDeterminantList(group1, group2, distance, iterative_interactions, versi def addIterativeAcidPair(object1, object2, interaction): - """ + """ Adding the Coulomb 'iterative' interaction (an acid pair): the higher pKa is raised with QQ+HB the lower pKa is lowered with HB - """ + """ values = interaction[1] annihilation = interaction[2] hbond_value = values[0] @@ -69,10 +74,10 @@ def addIterativeAcidPair(object1, object2, interaction): def addIterativeBasePair(object1, object2, interaction): - """ + """ Adding the Coulomb 'iterative' interaction (a base pair): the lower pKa is lowered - """ + """ values = interaction[1] annihilation = interaction[2] hbond_value = values[0] @@ -106,10 +111,10 @@ def addIterativeBasePair(object1, object2, interaction): def addIterativeIonPair(object1, object2, interaction, version): - """ + """ Adding the Coulomb 'iterative' interaction (an acid-base pair): the pKa of the acid is lowered & the pKa of the base is raised - """ + """ values = interaction[1] annihilation = interaction[2] hbond_value = values[0] @@ -132,7 +137,7 @@ def addIterativeIonPair(object1, object2, interaction, version): annihilation[0] = 0.00 annihilation[1] = 0.00 - + if add_term == True: # Coulomb @@ -161,9 +166,9 @@ def addIterativeIonPair(object1, object2, interaction, version): def addDeterminants(iterative_interactions, version, options=None): - """ + """ The iterative pKa scheme. Later it is all added in 'calculateTotalPKA' - """ + """ # --- setup --- iteratives = [] done_group = [] @@ -270,7 +275,7 @@ def addDeterminants(iterative_interactions, version, options=None): g = interaction[0] newDeterminant = Determinant(g, value) itres.group.determinants[type].append(newDeterminant) - + def findIterative(pair, iteratives): diff --git a/Source/lib.py b/propka/lib.py old mode 100755 new mode 100644 similarity index 91% rename from Source/lib.py rename to propka/lib.py index 43b471d..0bbb0a6 --- a/Source/lib.py +++ b/propka/lib.py @@ -1,23 +1,27 @@ -#!/usr/bin/python +from __future__ import division +from __future__ import print_function import string, sys, copy, math, os +import pkg_resources # # file I/O # def open_file_for_reading(filename): - if not os.path.isfile(filename): + try: + f = open(filename,'r') + except: raise Exception('Cannot find file %s' %filename) - - return open(filename,'r') + return f def open_file_for_writing(filename): - res = open(filename,'w') - if not res: + try: + res = open(filename,'w') + except: raise Exception('Could not open %s'%filename) return res - + # # bookkeeping etc. # @@ -43,7 +47,7 @@ def make_molecule(atom, atoms): if ba in atoms: atoms.remove(ba) res_atoms.extend(make_molecule(ba, atoms)) - + return res_atoms @@ -61,7 +65,7 @@ def generate_combinations(interactions): res.remove([]) return res - + def make_combination(combis, interaction): res = [] @@ -85,37 +89,37 @@ def loadOptions(): parser = OptionParser(usage) # loading the parser - parser.add_option("-f", "--file", action="append", dest="filenames", + parser.add_option("-f", "--file", action="append", dest="filenames", help="read data from , i.e. is added to arguments") - parser.add_option("-r", "--reference", dest="reference", default="neutral", + parser.add_option("-r", "--reference", dest="reference", default="neutral", help="setting which reference to use for stability calculations [neutral/low-pH]") - parser.add_option("-c", "--chain", action="append", dest="chains", + parser.add_option("-c", "--chain", action="append", dest="chains", help="creating the protein with only a specified chain, note, chains without ID are labeled 'A' [all]") - parser.add_option("-t", "--thermophile", action="append", dest="thermophiles", + parser.add_option("-t", "--thermophile", action="append", dest="thermophiles", help="defining a thermophile filename; usually used in 'alignment-mutations'") - parser.add_option("-a", "--alignment", action="append", dest="alignment", + parser.add_option("-a", "--alignment", action="append", dest="alignment", help="alignment file connecting and [.pir]") - parser.add_option("-m", "--mutation", action="append", dest="mutations", + parser.add_option("-m", "--mutation", action="append", dest="mutations", help="specifying mutation labels which is used to modify according to, e.g. N25R/N181D") - parser.add_option("-v", "--version", dest="version_label", default="Jan15", + parser.add_option("-v", "--version", dest="version_label", default="Jan15", help="specifying the sub-version of propka [Jan15/Dec19]") - parser.add_option("-p", "--parameters",dest="parameters", default="propka.cfg", - help="set the parameter file") - parser.add_option("-z", "--verbose", dest="verbose", action="store_true", default=True, + parser.add_option("-p", "--parameters",dest="parameters", default=pkg_resources.resource_filename(__name__, "propka.cfg"), + help="set the parameter file [%default]") + parser.add_option("-z", "--verbose", dest="verbose", action="store_true", default=True, help="sleep during calculations") parser.add_option("-q", "--quiet", dest="verbose", action="store_false", help="sleep during calculations") - parser.add_option("-s", "--silent", dest="verbose", action="store_false", + parser.add_option("-s", "--silent", dest="verbose", action="store_false", help="not activated yet") - parser.add_option("--verbosity", dest="verbosity", action="store_const", + parser.add_option("--verbosity", dest="verbosity", action="store_const", help="level of printout - not activated yet") - parser.add_option("-o", "--pH", dest="pH", type="float", default=7.0, + parser.add_option("-o", "--pH", dest="pH", type="float", default=7.0, help="setting pH-value used in e.g. stability calculations [7.0]") parser.add_option("-w", "--window", dest="window", nargs=3, type="float", default=(0.0, 14.0, 1.0), help="setting the pH-window to show e.g. stability profiles [0.0, 14.0, 1.0]") parser.add_option("-g", "--grid", dest="grid", nargs=3, type="float", default=(0.0, 14.0, 0.1), help="setting the pH-grid to calculate e.g. stability related properties [0.0, 14.0, 0.1]") - parser.add_option("--mutator", dest="mutator", + parser.add_option("--mutator", dest="mutator", help="setting approach for mutating [alignment/scwrl/jackal]") parser.add_option("--mutator-option", dest="mutator_options", action="append", help="setting property for mutator [e.g. type=\"side-chain\"]") @@ -156,7 +160,7 @@ def makeTidyAtomLabel(name,element): """ Returns a 'tidier' atom label for printing the new pdbfile """ - + if len(name)>4:# if longer than 4, just truncate the name label=name[0:4] elif len(name)==4:# if lenght is 4, otherwise use the name as it is diff --git a/Source/ligand.py b/propka/ligand.py old mode 100755 new mode 100644 similarity index 94% rename from Source/ligand.py rename to propka/ligand.py index d745d51..caf2fd9 --- a/Source/ligand.py +++ b/propka/ligand.py @@ -1,7 +1,12 @@ #!/usr/bin/python -import sys, Source.calculations -from Source.vector_algebra import * +from __future__ import division +from __future__ import print_function + +import sys + +import propka.calculations +from propka.vector_algebra import * all_sybyl_types = [ @@ -157,7 +162,7 @@ def assign_sybyl_type(atom): bonded_elements[i]=atom.bonded_atoms[i].element - + # Aromatic carbon/nitrogen if aromatic: for ra in ring_atoms: @@ -187,7 +192,7 @@ def assign_sybyl_type(atom): nitrogens = atom.get_bonded_elements('N') oxygens = atom.get_bonded_elements('O') if len(nitrogens)==1 and len(oxygens)==1: - C = atom + C = atom N = nitrogens[0] O = oxygens[0] @@ -199,9 +204,9 @@ def assign_sybyl_type(atom): set_type(C,'C.2') set_type(O,'O.2') return - - if atom.element=='C': + + if atom.element=='C': # check for carboxyl if len(atom.bonded_atoms)==3 and list(bonded_elements.values()).count('O')==2: i1 = list(bonded_elements.values()).index('O') @@ -213,11 +218,11 @@ def assign_sybyl_type(atom): return - + # sp carbon if len(atom.bonded_atoms)<=2: for b in atom.bonded_atoms: - if Source.calculations.squared_distance(atom, b) < max_C_triple_bond_squared: + if propka.calculations.squared_distance(atom, b) < max_C_triple_bond_squared: set_type(atom,'C.1') set_type(b,b.element+'.1') if atom.sybyl_assigned: @@ -232,7 +237,7 @@ def assign_sybyl_type(atom): if len(b.bonded_atoms)<3 or is_planar(b): set_type(b,'N.pl3') return - + # sp3 carbon set_type(atom, 'C.3') return @@ -244,11 +249,11 @@ def assign_sybyl_type(atom): if is_planar(atom.bonded_atoms[0]): set_type(atom,'N.pl3') return - + if planar: set_type(atom,'N.pl3') return - + set_type(atom,'N.3') return @@ -262,7 +267,7 @@ def assign_sybyl_type(atom): the_carbon = atom.bonded_atoms[0] if len(the_carbon.bonded_atoms)==3 and the_carbon.count_bonded_elements('O')==2: [O1,O2] = the_carbon.get_bonded_elements('O') - + if len(O1.bonded_atoms)==1 and len(O2.bonded_atoms)==1: set_type(O1,'O.co2-') set_type(O2,'O.co2') @@ -270,7 +275,7 @@ def assign_sybyl_type(atom): return # check for X=O - if Source.calculations.squared_distance(atom, atom.bonded_atoms[0]) < max_C_double_bond_squared: + if propka.calculations.squared_distance(atom, atom.bonded_atoms[0]) < max_C_double_bond_squared: set_type(atom,'O.2') if atom.bonded_atoms[0].element=='C': set_type(atom.bonded_atoms[0],'C.2') @@ -287,7 +292,7 @@ def assign_sybyl_type(atom): set_type(atom.bonded_atoms[i2],'O.2') set_type(atom,'S.o2') return - + # check for SO4 if list(bonded_elements.values()).count('O')==4: no_O2 = 0 @@ -297,7 +302,7 @@ def assign_sybyl_type(atom): no_O2+=1 else: set_type(atom.bonded_atoms[i],'O.3') - + set_type(atom,'S.3') @@ -315,8 +320,8 @@ def assign_sybyl_type(atom): # # find oxygens only bonded to current phosphorus # bonded_oxygens_1 = [o for o in bonded_oxygens if len(o.get_bonded_heavy_atoms())==1] # # find the closest oxygen ... -# closest_oxygen = min(bonded_oxygens_1, -# key= lambda o:Source.calculations.squared_distance(atom,o)) +# closest_oxygen = min(bonded_oxygens_1, +# key= lambda o:propka.calculations.squared_distance(atom,o)) # # ... and set it to O.2 # set_type(closest_oxygen,'O.2') @@ -378,20 +383,20 @@ def are_atoms_planar(atoms): v1 = vector(atom1=atoms[0], atom2=atoms[1]) v2 = vector(atom1=atoms[0], atom2=atoms[2]) n = (v1**v2).rescale(1.0) - + margin = 0.20 for b in atoms[3:]: v = vector(atom1=atoms[0], atom2=b).rescale(1.0) #print(atoms[0],abs(v*n) ) if abs(v*n)>margin: return False - + return True def is_aromatic_ring(atoms): if len(atoms)<5: return False - + for i in range(len(atoms)): if not are_atoms_planar(atoms[i:]+atoms[:i]): return False diff --git a/Source/ligand_pka_values.py b/propka/ligand_pka_values.py similarity index 89% rename from Source/ligand_pka_values.py rename to propka/ligand_pka_values.py index 591a77c..fcf799f 100644 --- a/Source/ligand_pka_values.py +++ b/propka/ligand_pka_values.py @@ -1,6 +1,9 @@ #!/usr/bin/env python -import Source.molecular_container, Source.calculations, Source.calculations, Source.parameters, Source.pdb, Source.lib, os, subprocess, sys +from __future__ import division +from __future__ import print_function + +import propka.molecular_container, propka.calculations, propka.calculations, propka.parameters, propka.pdb, propka.lib, os, subprocess, sys class ligand_pka_values: def __init__(self, parameters): @@ -14,14 +17,14 @@ class ligand_pka_values: print(self.molconvert) return - - + + def find_in_path(self, program): path = os.environ.get('PATH').split(os.pathsep) l = [i for i in filter(lambda loc: os.access(loc, os.F_OK), map(lambda dir: os.path.join(dir, program),path))] - + if len(l) == 0: print('Error: Could not find %s. Please make sure that it is found in the path.'%program) sys.exit(-1) @@ -29,14 +32,14 @@ class ligand_pka_values: return l[0] def get_marvin_pkas_for_pdb_file(self, file, no_pkas=10, min_pH =-10, max_pH=20): - molecule = Source.molecular_container.Molecular_container(file) + molecule = propka.molecular_container.Molecular_container(file) self.get_marvin_pkas_for_molecular_container(molecule, no_pkas=no_pkas, min_pH =min_pH, max_pH=max_pH) return def get_marvin_pkas_for_molecular_container(self, molecule, no_pkas=10, min_pH =-10, max_pH=20): for name in molecule.conformation_names: filename = '%s_%s'%(molecule.name,name) - self.get_marvin_pkas_for_conformation_container(molecule.conformations[name], name=filename, reuse=molecule.options.reuse_ligand_mol2_file, + self.get_marvin_pkas_for_conformation_container(molecule.conformations[name], name=filename, reuse=molecule.options.reuse_ligand_mol2_file, no_pkas=no_pkas, min_pH =min_pH, max_pH=max_pH) return @@ -50,7 +53,7 @@ class ligand_pka_values: def get_marvin_pkas_for_atoms(self, atoms, name='temp', reuse=False, no_pkas=10, min_pH =-10, max_pH=20): # do one molecule at the time so we don't confuse marvin - molecules = Source.lib.split_atoms_into_molecules(atoms) + molecules = propka.lib.split_atoms_into_molecules(atoms) for i in range(len(molecules)): filename = '%s_%d.mol2'%(name, i+1) self.get_marvin_pkas_for_molecule(molecules[i], filename=filename, reuse=reuse, no_pkas=no_pkas, min_pH =min_pH, max_pH=max_pH) @@ -61,18 +64,18 @@ class ligand_pka_values: def get_marvin_pkas_for_molecule(self, atoms, filename='__tmp_ligand.mol2', reuse=False, no_pkas=10, min_pH =-10, max_pH=20): # print out structure unless we are using user-modified structure if not reuse: - Source.pdb.write_mol2_for_atoms(atoms, filename) + propka.pdb.write_mol2_for_atoms(atoms, filename) # check that we actually have a file to work with if not os.path.isfile(filename): print('Warning: Didn\'t find a user-modified file \'%s\' - generating one'%filename) - Source.pdb.write_mol2_for_atoms(atoms, filename) + propka.pdb.write_mol2_for_atoms(atoms, filename) # Marvin # calculate pKa values options = 'pka -a %d -b %d --min %f --max %f -d large'%(no_pkas, no_pkas, min_pH, max_pH) - (output,errors) = subprocess.Popen([self.cxcalc, filename]+options.split(), + (output,errors) = subprocess.Popen([self.cxcalc, filename]+options.split(), stdout=subprocess.PIPE, stderr=subprocess.PIPE).communicate() if len(errors)>0: @@ -86,7 +89,7 @@ class ligand_pka_values: # extract calculated pkas indices,pkas,types = self.extract_pkas(output) - + # store calculated pka values for i in range(len(indices)): atoms[indices[i]].marvin_pka = pkas[i] @@ -104,10 +107,10 @@ class ligand_pka_values: values = values.split('\t') # format values - types = [tags[i][0] for i in range(1,len(tags)-1) if len(values)>i and values[i] != ''] + types = [tags[i][0] for i in range(1,len(tags)-1) if len(values)>i and values[i] != ''] indices = [int(a)-1 for a in values[-1].split(',') if a !=''] values = [float(v.replace(',','.')) for v in values[1:-1] if v != ''] - + if len(indices) != len(values) != len(types): raise Exception('Lengths of atoms and pka values mismatch') diff --git a/Source/molecular_container.py b/propka/molecular_container.py old mode 100755 new mode 100644 similarity index 86% rename from Source/molecular_container.py rename to propka/molecular_container.py index 83812f7..1665e14 --- a/Source/molecular_container.py +++ b/propka/molecular_container.py @@ -1,15 +1,19 @@ #!/usr/bin/python # -# Molecular container for storing all contents of pdb files +# Molecular container for storing all contents of pdb files # # +from __future__ import division +from __future__ import print_function -import os, Source.pdb, sys, Source.version, Source.output, Source.conformation_container, Source.group, Source.lib +import os, sys + +import propka.pdb, propka.version, propka.output, propka.conformation_container, propka.group, propka.lib class Molecular_container: def __init__(self, input_file, options=None): # printing out header before parsing input - Source.output.printHeader() + propka.output.printHeader() # set up some values self.options = options @@ -21,11 +25,11 @@ class Molecular_container: # set the version if options: - parameters = Source.parameters.Parameters(self.options.parameters) + parameters = propka.parameters.Parameters(self.options.parameters) else: - parameters = Source.parameters.Parameters('propka.cfg') + parameters = propka.parameters.Parameters('propka.cfg') try: - exec('self.version = Source.version.%s(parameters)'%parameters.version) + exec('self.version = propka.version.%s(parameters)'%parameters.version) except: raise Exception('Error: Version %s does not exist'%parameters.version) @@ -33,7 +37,7 @@ class Molecular_container: if input_file_extension[0:4] == '.pdb': # input is a pdb file # read in atoms and top up containers to make sure that all atoms are present in all conformations - [self.conformations, self.conformation_names] = Source.pdb.read_pdb(input_file, self.version.parameters,self) + [self.conformations, self.conformation_names] = propka.pdb.read_pdb(input_file, self.version.parameters,self) if len(self.conformations)==0: print('Error: The pdb file does not seems to contain any molecular conformations') sys.exit(-1) @@ -41,7 +45,7 @@ class Molecular_container: self.top_up_conformations() # make a structure precheck - Source.pdb.protein_precheck(self.conformations, self.conformation_names) + propka.pdb.protein_precheck(self.conformations, self.conformation_names) # set up atom bonding and protonation self.version.setup_bonding_and_protonation(self) @@ -58,12 +62,12 @@ class Molecular_container: # write out the input file filename = self.file.replace(input_file_extension,'.propka_input') - Source.pdb.write_input(self, filename) + propka.pdb.write_input(self, filename) elif input_file_extension == '.propka_input': #input is a propka_input file - [self.conformations, self.conformation_names] = Source.pdb.read_input(input_file, self.version.parameters, self) + [self.conformations, self.conformation_names] = propka.pdb.read_input(input_file, self.version.parameters, self) # Extract groups - this merely sets up the groups found in the input file self.extract_groups() @@ -103,15 +107,15 @@ class Molecular_container: """ Identify the groups needed for pKa calculation """ for name in self.conformation_names: self.conformations[name].extract_groups() - + return def additional_setup_when_reading_input_file(self): for name in self.conformation_names: self.conformations[name].additional_setup_when_reading_input_file() - + return - + def calculate_pka(self): # calculate for each conformation @@ -120,21 +124,21 @@ class Molecular_container: # find non-covalently coupled groups self.find_non_covalently_coupled_groups() - + # find the average of the conformations self.average_of_conformations() # print out the conformation-average results - Source.output.printResult(self, 'AVR', self.version.parameters) + propka.output.printResult(self, 'AVR', self.version.parameters) return def average_of_conformations(self): # make a new configuration to hold the average values - avr_conformation = Source.conformation_container.Conformation_container(name='average', - parameters=self.conformations[self.conformation_names[0]].parameters, + avr_conformation = propka.conformation_container.Conformation_container(name='average', + parameters=self.conformations[self.conformation_names[0]].parameters, molecular_container=self) - + for group in self.conformations[self.conformation_names[0]].get_titratable_groups_and_cysteine_bridges(): # new group to hold average values avr_group = group.clone() @@ -161,13 +165,13 @@ class Molecular_container: def write_pka(self, filename=None, reference="neutral", direction="folding", options=None): #for name in self.conformation_names: - # Source.output.writePKA(self, self.version.parameters, filename='%s_3.1_%s.pka'%(self.name, name), - # conformation=name,reference=reference, + # propka.output.writePKA(self, self.version.parameters, filename='%s_3.1_%s.pka'%(self.name, name), + # conformation=name,reference=reference, # direction=direction, options=options) # write out the average conformation filename=os.path.join('%s.pka'%(self.name)) - + # if the display_coupled_residues option is true, # write the results out to an alternative pka file if self.options.display_coupled_residues: @@ -176,8 +180,8 @@ class Molecular_container: if hasattr(self.version.parameters, 'output_file_tag') and len(self.version.parameters.output_file_tag)>0: filename=os.path.join('%s_%s.pka'%(self.name,self.version.parameters.output_file_tag)) - Source.output.writePKA(self, self.version.parameters, filename=filename, - conformation='AVR',reference=reference, + propka.output.writePKA(self, self.version.parameters, filename=filename, + conformation='AVR',reference=reference, direction=direction, options=options) return @@ -185,7 +189,7 @@ class Molecular_container: def getFoldingProfile(self, conformation='AVR',reference="neutral", direction="folding", grid=[0., 14., 0.1], options=None): # calculate stability profile profile = [] - for ph in Source.lib.make_grid(*grid): + for ph in propka.lib.make_grid(*grid): ddg = self.conformations[conformation].calculate_folding_energy( pH=ph, reference=reference) #print(ph,ddg) profile.append([ph, ddg]) @@ -196,29 +200,29 @@ class Molecular_container: opt = min(opt, point, key=lambda v:v[1]) # find values within 80 % of optimum - range_80pct = [None, None] + range_80pct = [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)] # find stability range stability_range = [None, None] - stable_values = [p[0] for p in profile if p[1]< 0.0] - + 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)] - + return profile, opt, range_80pct, stability_range def getChargeProfile(self, conformation='AVR', grid=[0., 14., .1]): charge_profile = [] - for ph in Source.lib.make_grid(*grid): + for ph in propka.lib.make_grid(*grid): q_unfolded, q_folded = self.conformations[conformation].calculate_charge(self.version.parameters, pH=ph) charge_profile.append([ph, q_unfolded, q_folded]) return charge_profile - + def getPI(self, conformation='AVR', grid=[0., 14., 1], iteration=0): #print('staring',grid, iteration) # search @@ -240,7 +244,7 @@ class Molecular_container: return pi_folded_value, pi_unfolded_value - + if __name__ == '__main__': input_file = sys.argv[1] diff --git a/Source/output.py b/propka/output.py similarity index 97% rename from Source/output.py rename to propka/output.py index 16a97e1..1521ea2 100644 --- a/Source/output.py +++ b/propka/output.py @@ -1,5 +1,10 @@ + +from __future__ import division +from __future__ import print_function + import sys -import Source.lib + +import propka.lib def printHeader(): @@ -17,7 +22,7 @@ def writePDB(protein, file=None, filename=None, include_hydrogens=False, options """ Write the residue to the new pdbfile """ - + if file == None: # opening file if not given if filename == None: @@ -124,7 +129,7 @@ def getDeterminantSection(protein, conformation, parameters): # printing determinants for chain in protein.conformations[conformation].chains: for residue_type in parameters.write_out_order: - groups = [g for g in protein.conformations[conformation].groups if g.atom.chainID == chain] + groups = [g for g in protein.conformations[conformation].groups if g.atom.chainID == chain] for group in groups: if group.residue_type == residue_type: str += "%s" % ( group.getDeterminantString(parameters.remove_penalised_group) ) @@ -158,8 +163,8 @@ def getFoldingProfileSection(protein, conformation='AVR', direction="folding", r str += "\n" str += "Free energy of %9s (kcal/mol) as a function of pH (using %s reference)\n" % (direction, reference) - profile, [pH_opt, dG_opt], [dG_min, dG_max], [pH_min, pH_max] = protein.getFoldingProfile(conformation=conformation, - reference=reference, + profile, [pH_opt, dG_opt], [dG_min, dG_max], [pH_min, pH_max] = protein.getFoldingProfile(conformation=conformation, + reference=reference, direction=direction, grid=[0., 14., 0.1], options=options) if profile == None: str += "Could not determine folding profile\n" @@ -183,7 +188,7 @@ def getFoldingProfileSection(protein, conformation='AVR', direction="folding", r str += "Could not determine the pH-range where the free energy is negative\n\n" else: str += "The free energy is negative in the range %4.1lf - %4.1lf\n\n" % (pH_min, pH_max) - + return str @@ -209,7 +214,7 @@ def getChargeProfileSection(protein, conformation='AVR', options=None): str += "Could not determine the pI\n\n" else: str += "The pI is %5.2lf (folded) and %5.2lf (unfolded)\n" % (pI_pro, pI_mod) - + return str @@ -300,8 +305,8 @@ def getReferencesHeader(): str += " Journal of Chemical Theory and Computation, 7(2):525-537 (2011)\n" str += " \n" str += " Improved Treatment of Ligands and Coupling Effects in Empirical Calculation\n" - str += " and Rationalization of pKa Values\n" - str += " Chresten R. Sondergaard, Mats H.M. Olsson, Michal Rostkowski, and Jan H. Jensen\n" + str += " and Rationalization of pKa Values\n" + str += " Chresten R. Sondergaard, Mats H.M. Olsson, Michal Rostkowski, and Jan H. Jensen\n" str += " Journal of Chemical Theory and Computation, (2011)\n" str += " \n" str += "-------------------------------------------------------------------------------------------------------\n" @@ -357,13 +362,13 @@ def getTheLine(): # Interaction maps def make_interaction_map(name, list, interaction): - """ Print out an interaction map named 'name' of the groups in 'list' + """ Print out an interaction map named 'name' of the groups in 'list' based on the function 'interaction' """ - + # return an empty string, if the list is empty if len(list)==0: return '' - + # for long list, use condensed formatting if len(list)>10: res = 'Condensed form:\n' @@ -377,7 +382,7 @@ def make_interaction_map(name, list, interaction): res = '%s\n%12s'%(name,'') for g in list: res += '%9s | '%g.label - + # do the map for g1 in list: res += '\n%-12s'%(g1.label) @@ -386,6 +391,6 @@ def make_interaction_map(name, list, interaction): if interaction(g1, g2): tag = ' X ' res += '%10s| '%tag - + return res - + diff --git a/Source/parameters.py b/propka/parameters.py similarity index 97% rename from Source/parameters.py rename to propka/parameters.py index 80d50f1..edf0c2a 100644 --- a/Source/parameters.py +++ b/propka/parameters.py @@ -1,7 +1,12 @@ + +from __future__ import division +from __future__ import print_function + import math -import Source.lib as lib +import propka.lib as lib import sys, os +import pkg_resources # names and types of all key words in configuration file matrices = ['interaction_matrix'] @@ -45,16 +50,15 @@ class Parameters: #self.print_interaction_parameters_latex() #####self.print_interactions_latex() #sys.exit(0) - - + + return def read_parameters(self, file): # try to locate the parameters file try: - path = os.path.dirname(__file__) - ifile = os.path.join(path,'../'+file) + ifile = pkg_resources.resource_filename(__name__, file) input = lib.open_file_for_reading(ifile) except: input = lib.open_file_for_reading(file) @@ -94,7 +98,7 @@ class Parameters: elif len(words)==3 and words[0] in string_dictionaries: self.parse_to_string_dictionary(words) - + #print(words) return @@ -226,25 +230,25 @@ class Parameters: print(""" Titratable: CG ARG -C2N ARG +C2N ARG N30 N+/LYS N31 N+/LYS N32 N+/LYS -N33 N+/LYS +N33 N+/LYS NAR HIS OCO COO OP TYR/SER? -SH CYS +SH CYS Non-titratable: NP1 AMD? OH ROH O3 ? -CL -F -NAM -N1 -O2 +CL +F +NAM +N1 +O2 """) return @@ -283,7 +287,7 @@ O2 'N1' :[], 'O2' :[]} - + s = """ \\begin{longtable}{lllll} \\caption{Ligand interaction parameters. For interactions not listed, the default value of %s is applied.} @@ -296,11 +300,11 @@ Group1 & Group2 & Interaction & c1 &c2 \\\\ \\multicolumn{5}{l}{\\emph{continued from the previous page}}\\\\ \\toprule -Group1 & Group2 & Interaction & c1 &c2 \\\\ +Group1 & Group2 & Interaction & c1 &c2 \\\\ \\midrule \\endhead -\\midrule +\\midrule \\multicolumn{5}{r}{\\emph{continued on the next page}}\\\\ \\endfoot @@ -332,7 +336,7 @@ Group1 & Group2 & Interaction & c1 &c2 \\\\ agroups = ['COO', 'HIS', 'CYS', 'TYR', 'SER', 'N+', 'LYS', 'AMD', 'ARG', 'TRP', 'ROH', 'CG', 'C2N', 'N30', 'N31', 'N32', 'N33', 'NAR', 'OCO', 'NP1', 'OH', 'O3', 'CL', 'F', 'NAM', 'N1', 'O2', 'OP', 'SH'] lgroups = ['CG', 'C2N', 'N30', 'N31', 'N32', 'N33', 'NAR', 'OCO', 'NP1', 'OH', 'O3', 'CL', 'F', 'NAM', 'N1', 'O2', 'OP', 'SH'] - + s = """ \\begin{longtable}{%s} \\caption{Ligand interaction parameters. For interactions not listed, the default value of %s is applied.} @@ -345,11 +349,11 @@ Group1 & Group2 & Interaction & c1 &c2 \\\\ \\multicolumn{5}{l}{\\emph{continued from the previous page}}\\\\ \\toprule -Group1 & Group2 & Interaction & c1 &c2 \\\\ +Group1 & Group2 & Interaction & c1 &c2 \\\\ \\midrule \\endhead -\\midrule +\\midrule \\multicolumn{5}{r}{\\emph{continued on the next page}}\\\\ \\endfoot @@ -399,14 +403,14 @@ class Interaction_matrix: self.dictionary[group][new_group] = self.value self.dictionary[new_group][group] = self.value - + return def get_value(self, item1, item2): try: return self.dictionary[item1][item2] except: - return None + return None def __getitem__(self, group): if group not in self.dictionary.keys(): @@ -430,7 +434,7 @@ class Interaction_matrix: return s # ks = ['COO', 'SER', 'ARG', 'LYS', 'HIS', 'AMD', 'CYS', 'TRP','ROH','TYR','N+','CG', 'C2N', 'N30', 'N31', 'N32', 'N33', 'NAR', 'OCO', 'NP1', 'OH', 'O3', 'CL', 'F', 'NAM', 'N1', 'O2', 'OP', 'SH'] - + # p = '' # n=0 # for i in range(len(ks)): @@ -442,7 +446,7 @@ class Interaction_matrix: # print('total',n,len(ks)) # return p - + class Pair_wise_matrix: @@ -451,7 +455,7 @@ class Pair_wise_matrix: self.dictionary = {} self.default = [0.0, 0.0] return - + def add(self,words): # assign the default value if len(words)==3 and words[0]=='default': @@ -462,7 +466,7 @@ class Pair_wise_matrix: g1 = words[0] g2 = words[1] v = [float(words[2]), float(words[3])] - + self.insert(g1,g2,v) self.insert(g2,g1,v) @@ -482,7 +486,7 @@ class Pair_wise_matrix: return def get_value(self, item1, item2): - + try: return self.dictionary[item1][item2] except: @@ -504,7 +508,7 @@ class Pair_wise_matrix: s += '%s %s %s\n'%(k1,k2,self[k1][k2]) return s - + diff --git a/Source/pdb.py b/propka/pdb.py similarity index 91% rename from Source/pdb.py rename to propka/pdb.py index d0a3bfa..9443b35 100644 --- a/Source/pdb.py +++ b/propka/pdb.py @@ -1,6 +1,12 @@ -import string, sys, copy, Source.lib -from Source.atom import Atom -from Source.conformation_container import Conformation_container + +from __future__ import division +from __future__ import print_function + +import string, sys, copy + +import propka.lib +from propka.atom import Atom +from propka.conformation_container import Conformation_container expected_atom_numbers = {'ALA':5, 'ARG':11, @@ -35,23 +41,23 @@ def read_pdb(pdb_file, parameters, molecule): conformations[name].add_atom(atom) # make a sorted list of conformation names - names = sorted(conformations.keys(), key=Source.lib.conformation_sorter) - + names = sorted(conformations.keys(), key=propka.lib.conformation_sorter) + return [conformations, names] def protein_precheck(conformations, names): - + for name in names: atoms = conformations[name].atoms - + res_ids = [] [res_ids.append(resid_from_atom(a)) for a in atoms if not res_ids.count(resid_from_atom(a))] - + for res_id in res_ids: res_atoms = [a for a in atoms if resid_from_atom(a) == res_id and a.element != 'H'] resname = res_atoms[0].resName residue_label = '%3s%5s'%(resname, res_id) - + # ignore ligand residues if resname not in expected_atom_numbers: continue @@ -65,7 +71,7 @@ def protein_precheck(conformations, names): # check number of atoms in residue if len(res_atoms) != expected_atom_numbers[resname]: print('Warning: Unexpected number (%d) of atoms in residue %s in conformation %s'%(len(res_atoms),residue_label, name)) - + return def resid_from_atom(a): @@ -74,7 +80,7 @@ def resid_from_atom(a): def get_atom_lines_from_pdb(pdb_file, ignore_residues = [], keep_protons=False, tags = ['ATOM ', 'HETATM'], chains=None): - lines = Source.lib.open_file_for_reading(pdb_file).readlines() + lines = propka.lib.open_file_for_reading(pdb_file).readlines() nterm_residue = 'next_residue' old_residue = None terminal = None @@ -88,7 +94,7 @@ def get_atom_lines_from_pdb(pdb_file, ignore_residues = [], keep_protons=False, if tag == 'MODEL ': model = int(line[6:]) nterm_residue = 'next_residue' - + if tag == 'TER ': nterm_residue = 'next_residue' @@ -102,21 +108,21 @@ def get_atom_lines_from_pdb(pdb_file, ignore_residues = [], keep_protons=False, continue if chains and line[21] not in chains: continue - - # set the Nterm residue number - nessecary because we may need to + + # set the Nterm residue number - nessecary because we may need to # identify more than one N+ group for structures with alt_conf tags if nterm_residue == 'next_residue' and tag == 'ATOM ': - # make sure that we reached a new residue - nessecary if OXT is not the last atom in - # the previous residue + # make sure that we reached a new residue - nessecary if OXT is not the last atom in + # the previous residue if old_residue != residue_number: nterm_residue = residue_number old_residue = None - # Identify the configuration + # Identify the configuration # convert digits to letters if alt_conf_tag in '123456789': - alt_conf_tag = chr(ord(alt_conf_tag)+16) + alt_conf_tag = chr(ord(alt_conf_tag)+16) if alt_conf_tag == ' ': alt_conf_tag = 'A' conformation = '%d%s'%(model, alt_conf_tag) @@ -145,10 +151,10 @@ def get_atom_lines_from_pdb(pdb_file, ignore_residues = [], keep_protons=False, def write_pdb(conformation, filename): write_pdb_for_atoms(conformation.atoms, filename) return - + def write_pdb_for_atoms(atoms, filename, make_conect_section=False): - out = Source.lib.open_file_for_writing(filename) - + out = propka.lib.open_file_for_writing(filename) + for atom in atoms: out.write(atom.make_pdb_line()) @@ -171,7 +177,7 @@ def write_mol2_for_atoms(atoms, filename): for i in range(len(atoms)): atoms_section += atoms[i].make_mol2_line(i+1) - + bonds_section = '@BOND\n' id = 1 for i in range(len(atoms)): @@ -180,12 +186,12 @@ def write_mol2_for_atoms(atoms, filename): type = get_bond_order(atoms[i],atoms[j]) bonds_section += '%7d %7d %7d %7s\n'%(id, i+1, j+1, type) id+=1 - + substructure_section = '@SUBSTRUCTURE\n\n' if len(atoms)>0: substructure_section = '@SUBSTRUCTURE\n%-7d %10s %7d\n'%(atoms[0].resNumb,atoms[0].resName,atoms[0].numb) - out = Source.lib.open_file_for_writing(filename) + out = propka.lib.open_file_for_writing(filename) out.write(header%(len(atoms),id-1)) out.write(atoms_section) out.write(bonds_section) @@ -209,14 +215,14 @@ def get_bond_order(atom1, atom2): if '.ar' in atom1.sybyl_type and '.ar' in atom2.sybyl_type: type = 'ar' - + return type def write_input(molecular_container, filename): - out = Source.lib.open_file_for_writing(filename) + out = propka.lib.open_file_for_writing(filename) for conformation_name in molecular_container.conformation_names: out.write('MODEL %s\n'%conformation_name) @@ -250,14 +256,14 @@ def read_input(input_file, parameters,molecule): conformations[name].add_atom(atom) # make a sorted list of conformation names - names = sorted(conformations.keys(), key=Source.lib.conformation_sorter) - + names = sorted(conformations.keys(), key=propka.lib.conformation_sorter) + return [conformations, names] def get_atom_lines_from_input(input_file, tags = ['ATOM ','HETATM']): - lines = Source.lib.open_file_for_reading(input_file).readlines() + lines = propka.lib.open_file_for_reading(input_file).readlines() conformation = '' atoms = {} @@ -278,7 +284,7 @@ def get_atom_lines_from_input(input_file, tags = ['ATOM ','HETATM']): atom.is_protonated = True atoms[atom.numb] = atom numbers.append(atom.numb) - + # found bonding information - apply it if tag == 'CONECT' and len(line)>14: conect_numbers = [line[i:i+5] for i in range(6, len(line)-1, 5)] @@ -301,7 +307,7 @@ def get_atom_lines_from_input(input_file, tags = ['ATOM ','HETATM']): center_atom = atoms[int(conect_numbers[0])] for n in conect_numbers[1:]: cg = atoms[int(n)] - center_atom.group.couple_covalently(cg.group) + center_atom.group.couple_covalently(cg.group) # found info on non-covalent coupling if tag == 'NCOUPL' and len(line)>14: @@ -309,8 +315,8 @@ def get_atom_lines_from_input(input_file, tags = ['ATOM ','HETATM']): center_atom = atoms[int(conect_numbers[0])] for n in conect_numbers[1:]: cg = atoms[int(n)] - center_atom.group.couple_non_covalently(cg.group) - + center_atom.group.couple_non_covalently(cg.group) + # this conformation is done - yield the atoms if tag == 'ENDMDL': for n in numbers: diff --git a/propka/propka.cfg b/propka/propka.cfg new file mode 100644 index 0000000..a12e38a --- /dev/null +++ b/propka/propka.cfg @@ -0,0 +1,398 @@ +# PropKa configuration file + +version version_A + +# Model pKa values +model_pkas C- 3.20 +model_pkas ASP 3.80 +model_pkas GLU 4.50 +model_pkas HIS 6.50 +model_pkas CYS 9.00 +model_pkas TYR 10.00 +model_pkas LYS 10.50 +model_pkas ARG 12.50 +#model_pkas SER 14.20 Jack Kyte: Structure in Protein Chemistry, 1995, Garland Publishing, Inc New York and London +model_pkas N+ 8.00 +model_pkas CG 11.50 +model_pkas C2N 11.50 +model_pkas N30 10.00 +model_pkas N31 10.00 +model_pkas N32 10.00 +model_pkas N33 10.00 +model_pkas NAR 5.00 +model_pkas OCO 4.50 +model_pkas SH 10.00 +model_pkas OP 6.00 + +# Custom ligand pKa values +# P. Acharya, P. Cheruku, S. Chatterjee, S. Acharya, and, J. Chattopadhyaya: +# Measurement of Nucleobase pKa Values in Model Mononucleotides +# Shows RNA-RNA Duplexes To Be More Stable than DNA-DNA Duplexes +# Journal of the American Chemical Society 2004 126 (9), 2862-2869 +# +custom_model_pkas DA-N1 3.82 +custom_model_pkas DA-N3 3.82 +custom_model_pkas DA-N7 3.82 +custom_model_pkas DA-OP1 1.00 +custom_model_pkas DA-OP2 1.00 + +custom_model_pkas DG-N1 9.59 +custom_model_pkas DG-N3 9.59 +custom_model_pkas DG-N7 9.59 +custom_model_pkas DG-OP1 1.00 +custom_model_pkas DG-OP2 1.00 + +custom_model_pkas DC-N3 4.34 +custom_model_pkas DC-OP1 1.00 +custom_model_pkas DC-OP2 1.00 + +custom_model_pkas DT-N3 10.12 +custom_model_pkas DT-OP1 1.00 +custom_model_pkas DT-OP2 1.00 + + +# protein group mapping +protein_group_mapping ASP-CG COO +protein_group_mapping GLU-CD COO +protein_group_mapping HIS-CG HIS +protein_group_mapping CYS-SG CYS +protein_group_mapping TYR-OH TYR +protein_group_mapping LYS-NZ LYS +protein_group_mapping ARG-CZ ARG +#protein_group_mapping SER-OG SER +protein_group_mapping THR-OG1 ROH +protein_group_mapping SER-OG ROH# +protein_group_mapping ASN-CG AMD +protein_group_mapping GLN-CD AMD +protein_group_mapping TRP-NE1 TRP + + +# matrix for propka interactions +# 'N' non-iterative interaction +# 'I' iterative interaction +# '-' no interaction + #CYS +interaction_matrix CYS I#N+ +interaction_matrix N+ N I#HIS +interaction_matrix HIS I N I#LYS +interaction_matrix LYS N N N I#AMD +interaction_matrix AMD N - N - -#COO +interaction_matrix COO I N I N N I#ARG +interaction_matrix ARG N N N N - N I#TRP +interaction_matrix TRP N - - - - N - -#ROH +interaction_matrix ROH N - - - - N - - -#TYR +interaction_matrix TYR N I I I N N N N N I#SER +interaction_matrix SER N N N N N N I N N N I #CG +interaction_matrix CG N N N N - N I - - N I I#C2N +interaction_matrix C2N N N N N - N I - - N I I I#N30 +interaction_matrix N30 N I N N - N N - - I N I I I#N31 +interaction_matrix N31 N I N N - N N - - I N I I I I#N32 +interaction_matrix N32 N I N N - N N - - I N I I I I I#N33 +interaction_matrix N33 N I N N - N N - - I N I I I I I I#NAR +interaction_matrix NAR I N I I N I N - - I N N N N N N N I#OCO +interaction_matrix OCO I N I N N I N N N N N N N N N N N I I#NP1 +interaction_matrix NP1 N - N - - N - - - N N - - - - - - N N -#OH +interaction_matrix OH N - - - - N - - - N N - - - - - - - N - -#O3 +interaction_matrix O3 N - N - - N - - - N N - - - - - - N N - - -#CL +interaction_matrix CL N - N - - N - - - N N - - - - - - N N - - - -#F +interaction_matrix F N - N - - N - - - N N - - - - - - N N - - - - -#NAM +interaction_matrix NAM N - N - - N - - - N N - - - - - - N N - - - - - -#N1 +interaction_matrix N1 N - N - - N - - - N N - - - - - - N N - - - - - - -#O2 +interaction_matrix O2 N - N - - N - - - N N - - - - - - N N - - - - - - - -#OP +interaction_matrix OP I N I N N I N N N N N N N N N N N I I N N N N N N N N I#SH +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 + +# Cutoff values for side chain interactions +# default value +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 ARG 1.85 2.85 +sidechain_cutoffs COO LYS 2.85 3.85 +sidechain_cutoffs COO HIS 2.0 3.0 +sidechain_cutoffs COO AMD 2.0 3.0 +sidechain_cutoffs COO TRP 2.0 3.0 +sidechain_cutoffs COO ROH 2.65 3.65 +sidechain_cutoffs COO TYR 2.65 3.65 +sidechain_cutoffs COO N+ 2.85 3.85 +sidechain_cutoffs COO CG 1.85 2.85 +sidechain_cutoffs COO C2N 1.85 2.85 +sidechain_cutoffs COO N30 2.85 3.85 +sidechain_cutoffs COO N31 2.85 3.85 +sidechain_cutoffs COO N32 2.85 3.85 +sidechain_cutoffs COO N33 2.85 3.85 +sidechain_cutoffs COO NAR 2.0 3.0 +sidechain_cutoffs COO OCO 2.5 3.5 +sidechain_cutoffs COO OH 2.65 3.65 +sidechain_cutoffs COO NAM 2.0 3.0 +# SER +sidechain_cutoffs SER SER 3.5 4.5 +sidechain_cutoffs SER ARG 2.5 4.0 +sidechain_cutoffs SER HIS 2.0 3.0 +sidechain_cutoffs SER AMD 2.5 3.5 +sidechain_cutoffs SER CYS 3.5 4.5 +sidechain_cutoffs SER TRP 2.5 3.5 +sidechain_cutoffs SER ROH 3.5 4.5 +sidechain_cutoffs SER CG 2.5 4.0 +sidechain_cutoffs SER C2N 2.5 4.0 +sidechain_cutoffs SER NAR 2.0 3.0 +sidechain_cutoffs SER OH 3.5 4.5 +sidechain_cutoffs SER SH 3.5 4.5 +sidechain_cutoffs SER TYR 3.5 4.5 +sidechain_cutoffs SER N+ 3.0 4.5 +sidechain_cutoffs SER NAM 2.5 3.5 +# ARG +sidechain_cutoffs ARG CYS 2.5 4.0 +sidechain_cutoffs ARG TYR 2.5 4.0 +sidechain_cutoffs ARG OCO 1.85 2.85 +sidechain_cutoffs ARG SH 2.5 4.0 +# HIS +sidechain_cutoffs HIS AMD 2.0 3.0 +sidechain_cutoffs HIS TYR 2.0 3.0 +sidechain_cutoffs HIS OCO 2.0 3.0 +# CYS +sidechain_cutoffs CYS CYS 3.0 5.0 +sidechain_cutoffs CYS TRP 2.5 3.5 +sidechain_cutoffs CYS ROH 3.5 4.5 +sidechain_cutoffs CYS AMD 2.5 3.5 +sidechain_cutoffs CYS TYR 3.5 4.5 +sidechain_cutoffs CYS N+ 3.0 4.5 +sidechain_cutoffs CYS CG 2.5 4.0 +sidechain_cutoffs CYS C2N 2.5 4.0 +sidechain_cutoffs CYS N30 3.0 4.5 +sidechain_cutoffs CYS N31 3.0 4.5 +sidechain_cutoffs CYS N32 3.0 4.5 +sidechain_cutoffs CYS N33 3.0 4.5 +sidechain_cutoffs CYS OH 3.5 4.5 +sidechain_cutoffs CYS NAM 2.5 3.5 +sidechain_cutoffs CYS SH 3.0 5.0 +# TYR +sidechain_cutoffs TYR TYR 3.5 4.5 +sidechain_cutoffs TYR N+ 3.0 4.5 +sidechain_cutoffs TYR AMD 2.5 3.5 +sidechain_cutoffs TYR TRP 2.5 3.5 +sidechain_cutoffs TYR ROH 3.5 4.5 +sidechain_cutoffs TYR CG 2.5 4.0 +sidechain_cutoffs TYR C2N 2.5 4.0 +sidechain_cutoffs TYR OCO 2.65 3.65 +sidechain_cutoffs TYR NAR 2.0 3.0 +sidechain_cutoffs TYR OH 3.5 4.5 +sidechain_cutoffs TYR NAM 2.5 3.5 +sidechain_cutoffs TYR SH 3.5 4.5 +# N+ +sidechain_cutoffs N+ OCO 2.85 3.85 +sidechain_cutoffs N+ SH 3.0 4.5 +# LYS +sidechain_cutoffs LYS OCO 2.85 3.85 +# OCO +sidechain_cutoffs OCO OCO 2.5 3.5 +sidechain_cutoffs OCO TRP 2.0 3.0 +sidechain_cutoffs OCO ROH 2.65 3.65 +sidechain_cutoffs OCO AMD 2.0 3.0 +sidechain_cutoffs OCO CG 1.85 2.85 +sidechain_cutoffs OCO C2N 1.85 2.85 +sidechain_cutoffs OCO N30 2.85 3.85 +sidechain_cutoffs OCO N31 2.85 3.85 +sidechain_cutoffs OCO N32 2.85 3.85 +sidechain_cutoffs OCO N33 2.85 3.85 +sidechain_cutoffs OCO NAR 2.0 3.0 +sidechain_cutoffs OCO OH 2.65 3.65 +sidechain_cutoffs OCO NAM 2.0 3.0 +# NAR +sidechain_cutoffs NAR AMD 2.0 3.0 +# SH +sidechain_cutoffs SH ROH 3.5 4.5 +sidechain_cutoffs SH TRP 2.5 3.5 +sidechain_cutoffs SH AMD 2.5 3.5 +sidechain_cutoffs SH NAM 2.5 3.5 +sidechain_cutoffs SH CG 2.5 4.0 +sidechain_cutoffs SH C2N 2.5 4.0 +sidechain_cutoffs SH OH 3.5 4.5 +sidechain_cutoffs SH SH 3.0 5.0 + + + +# Maximal interaction energies for side chains +sidechain_interaction 0.85 + +# Angular dependent sidechain interactions +angular_dependent_sidechain_interactions HIS +angular_dependent_sidechain_interactions ARG +angular_dependent_sidechain_interactions AMD +angular_dependent_sidechain_interactions TRP + +# exception interaction values +COO_HIS_exception 1.60 +OCO_HIS_exception 1.60 +CYS_HIS_exception 1.60 +CYS_CYS_exception 3.60 + +# Coulomb interaction parameters +coulomb_cutoff1 4.0 +coulomb_cutoff2 10.0 +coulomb_diel 80.0 + +# Backbone hydrogen bond parameters +backbone_NH_hydrogen_bond COO -0.85 2.00 3.00 +#backbone_NH_hydrogen_bond C- -0.85 2.00 3.00 +backbone_NH_hydrogen_bond CYS -0.85 3.00 4.00 +backbone_NH_hydrogen_bond TYR -0.85 2.20 3.20 +backbone_NH_hydrogen_bond OCO -0.85 2.00 3.50 +backbone_NH_hydrogen_bond NAR -0.85 2.00 3.50 + +backbone_CO_hydrogen_bond HIS 0.85 2.00 3.00 +backbone_CO_hydrogen_bond OCO 0.85 3.00 4.00 +backbone_CO_hydrogen_bond CG 0.85 2.00 4.00 +backbone_CO_hydrogen_bond C2N 0.85 2.00 4.00 +backbone_CO_hydrogen_bond N30 0.85 2.00 4.00 +backbone_CO_hydrogen_bond N31 0.85 2.00 4.00 +backbone_CO_hydrogen_bond N32 0.85 2.00 4.00 +backbone_CO_hydrogen_bond N33 0.85 2.00 4.00 +backbone_CO_hydrogen_bond NAR 0.85 2.00 3.50 + +# Group charges +charge COO -1 +charge HIS +1 +charge CYS -1 +charge TYR -1 +charge LYS +1 +charge ARG +1 +charge N+ +1 +charge C- -1 +charge OCO -1 +charge SER -1 +charge CG +1 +charge C2N +1 +charge N30 +1 +charge N31 +1 +charge N32 +1 +charge N33 +1 +charge NAR +1 +charge SH -1 +charge OP -1 + +# list of acids +acid_list ASP +acid_list GLU +acid_list CYS +acid_list TYR +acid_list SER +acid_list C- +acid_list OCO +acid_list OP +acid_list SH + +# list of bases +base_list ARG +base_list LYS +base_list HIS +base_list N+ +base_list CG +base_list C2N +base_list N30 +base_list N31 +base_list N32 +base_list N33 +base_list NAR + +# list of groups used in backbone reorganisation calculations +backbone_reorganisation_list ASP +backbone_reorganisation_list GLU + +# Residues that should be ignored +ignore_residues HOH +ignore_residues H2O +ignore_residues HOH +ignore_residues SO4 +ignore_residues PO4 +ignore_residues PEG +ignore_residues EPE +#ignore_residues NAG +ignore_residues TRS + +# Relative Van der Waals volume parameters for the radial volume model +# Radii adopted from Bondi, A. (1964). "Van der Waals Volumes and Radii". J. Phys. Chem. 68 (3): 441-51 +VanDerWaalsVolume C 1.40 # radius: 1.70, volume: 20.58 all 'C' and 'CA' atoms +VanDerWaalsVolume C4 2.64 # 38.79 hydrodphobic carbon atoms + unidentified atoms +VanDerWaalsVolume N 1.06 # radius: 1.55, volume: 15.60 all nitrogen atoms +VanDerWaalsVolume O 1.00 # radius: 1.52, volume: 14.71 all oxygen atoms +VanDerWaalsVolume S 1.66 # radius: 1.80, volume: 24.43 all sulphur atoms +VanDerWaalsVolume F 0.90 # raidus: 1.47, volume: 13.30 for fluorine +VanDerWaalsVolume Cl 1.53 # radius: 1.75, volume: 22.44 for chlorine +VanDerWaalsVolume P 1.66 # radius: 1.80, volume: 24.42 for phosphorus + +# Other desolvation parameters +desolvationSurfaceScalingFactor 0.25 +desolvationPrefactor -13.0 +desolvationAllowance 0.0 +desolv_cutoff 20.0 +buried_cutoff 15.0 +Nmin 280 +Nmax 560 + +# Ligand groups +ligand_typing groups +min_bond_distance_for_hydrogen_bonds 4 + +# covalent coupling +coupling_max_number_of_bonds 3 +shared_determinants 0 +common_charge_centre 0 +remove_penalised_group 1 + +# non-covalent coupling +max_intrinsic_pKa_diff 2.0 +min_interaction_energy 0.5 +max_free_energy_diff 1.0 +min_swap_pka_shift 1.0 +min_pka 0.0 +max_pka 10.0 +pH variable +reference neutral + +# ions +ions 1P 1 # generic charged atoms +ions 2P 2 +ions 1N -1 +ions 2N -2 + +ions MG 2 #Magnesium Ion +ions CA 2 #Calcium Ion +ions ZN 2 #Zinc Ion +ions NA 1 #Sodium Ion +ions CL -1 #Chloride Ion +ions MN 2 #Manganese (ii) Ion +ions K 1 #Potassium Ion +ions CD 2 #Cadmium Ion +ions FE 3 #Fe (iii) Ion +ions SR 2 #Strontium Ion +ions CU 2 #Copper (ii) Ion +ions IOD -1 #Iodide Ion +ions HG 2 #Mercury (ii) Ion +ions BR -1 #Bromide Ion +ions CO 2 #Cobalt (ii) Ion +ions NI 2 #Nickel (ii) Ion +ions FE2 2 #Fe (ii) Ion + +# write out order of residues +write_out_order ASP +write_out_order GLU +write_out_order C- +write_out_order HIS +write_out_order CYS +write_out_order TYR +write_out_order LYS +write_out_order ARG +write_out_order SER +write_out_order N+ +write_out_order CG +write_out_order C2N +write_out_order N30 +write_out_order N31 +write_out_order N32 +write_out_order N33 +write_out_order NAR +write_out_order OCO +write_out_order SH +write_out_order OP diff --git a/Source/protein_bonds.dat b/propka/protein_bonds.dat similarity index 100% rename from Source/protein_bonds.dat rename to propka/protein_bonds.dat diff --git a/Source/protonate.py b/propka/protonate.py old mode 100755 new mode 100644 similarity index 93% rename from Source/protonate.py rename to propka/protonate.py index 92ab0f9..396790b --- a/Source/protonate.py +++ b/propka/protonate.py @@ -1,11 +1,14 @@ #!/usr/bin/python -from Source.vector_algebra import * -import Source.bonds, Source.pdb, Source.atom +from __future__ import division +from __future__ import print_function + +from propka.vector_algebra import * +import propka.bonds, propka.pdb, propka.atom class Protonate: """ Protonates atoms using VSEPR theory """ - + def __init__(self, verbose=False): self.verbose=verbose @@ -56,9 +59,9 @@ class Protonate: 'HIS-ND1':1.0, 'LYS-NZ':1.0, 'N+':1.0, - 'C-':-1.0} + 'C-':-1.0} + - self.sybyl_charges = {'N.pl3':+1, 'N.3':+1, 'N.4':+1, @@ -92,17 +95,17 @@ class Protonate: self.display('----- Protonation started -----') # Remove all currently present hydrogen atoms self.remove_all_hydrogen_atoms(molecules) - + # protonate all atoms for name in molecules.conformation_names: non_H_atoms = molecules.conformations[name].get_non_hydrogen_atoms() for atom in non_H_atoms: self.protonate_atom(atom) - + # fix hydrogen names #self.set_proton_names(non_H_atoms) - + return @@ -110,7 +113,7 @@ class Protonate: for name in molecular_container.conformation_names: molecular_container.conformations[name].atoms = molecular_container.conformations[name].get_non_hydrogen_atoms() return - + def set_charge(self, atom): # atom is a protein atom @@ -148,11 +151,11 @@ class Protonate: for heavy_atom in heavy_atoms: i = 1 for bonded in heavy_atom.bonded_atoms: - + if bonded.element == 'H': bonded.name+='%d'%i i+=1 - + return @@ -160,7 +163,7 @@ class Protonate: def set_number_of_protons_to_add(self, atom): self.display('*'*10) self.display('Setting number of protons to add for',atom) - atom.number_of_protons_to_add = 8 + atom.number_of_protons_to_add = 8 self.display(' %4d'%8) atom.number_of_protons_to_add -= self.valence_electrons[atom.element] self.display('Valence eletrons: %4d'%-self.valence_electrons[atom.element]) @@ -177,7 +180,7 @@ class Protonate: return def set_steric_number_and_lone_pairs(self, atom): - + # If we already did this, there is no reason to do it again if atom.steric_number_and_lone_pairs_set: return @@ -195,10 +198,10 @@ class Protonate: atom.steric_number = 0 - + self.display('%65s: %4d'%('Valence electrons',self.valence_electrons[atom.element])) atom.steric_number += self.valence_electrons[atom.element] - + self.display('%65s: %4d'%('Number of bonds',len(atom.bonded_atoms))) atom.steric_number += len(atom.bonded_atoms) @@ -216,7 +219,7 @@ class Protonate: self.display('%65s: %4.1f'%('Charge(-)',atom.charge)) atom.steric_number -= atom.charge - + atom.steric_number = math.floor(atom.steric_number/2.0) atom.number_of_lone_pairs = atom.steric_number - len(atom.bonded_atoms) - atom.number_of_protons_to_add @@ -240,9 +243,9 @@ class Protonate: return - + def trigonal(self, atom): - self.display('TRIGONAL - %d bonded atoms'%(len(atom.bonded_atoms))) + self.display('TRIGONAL - %d bonded atoms'%(len(atom.bonded_atoms))) rot_angle = math.radians(120.0) c = vector(atom1 = atom) @@ -250,13 +253,13 @@ class Protonate: # 0 bonds if len(atom.bonded_atoms) == 0: pass - + # 1 bond if len(atom.bonded_atoms) == 1 and atom.number_of_protons_to_add > 0: # Add another atom with the right angle to the first one a = vector(atom1 = atom, atom2 = atom.bonded_atoms[0]) # use plane of bonded trigonal atom - e.g. arg - + self.set_steric_number_and_lone_pairs(atom.bonded_atoms[0]) if atom.bonded_atoms[0].steric_number == 3 and len(atom.bonded_atoms[0].bonded_atoms)>1: # use other atoms bonded to the neighbour to establish the plane, if possible @@ -265,21 +268,21 @@ class Protonate: if atom.bonded_atoms[0].bonded_atoms[i] != atom: other_atom_indices.append(i) - + v1 = vector(atom1 = atom, atom2 = atom.bonded_atoms[0]) - v2 = vector(atom1 = atom.bonded_atoms[0], + v2 = vector(atom1 = atom.bonded_atoms[0], atom2 = atom.bonded_atoms[0].bonded_atoms[other_atom_indices[0]]) - + axis = v1**v2 - + # 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: - v3 = vector(atom1 = atom.bonded_atoms[0], + v3 = vector(atom1 = atom.bonded_atoms[0], atom2 = atom.bonded_atoms[0].bonded_atoms[other_atom_indices[1]]) axis2 = v1**v3 - + if axis * axis2>0: axis = axis+axis2 else: @@ -294,12 +297,12 @@ class Protonate: # 2 bonds if len(atom.bonded_atoms) == 2 and atom.number_of_protons_to_add > 0: - # Add another atom with the right angle to the first two + # Add another atom with the right angle to the first two a1 = vector(atom1 = atom, atom2 = atom.bonded_atoms[0]).rescale(1.0) a2 = vector(atom1 = atom, atom2 = atom.bonded_atoms[1]).rescale(1.0) - - new_a = -a1 - a2 - new_a = self.set_bond_distance(new_a, atom.element) + + new_a = -a1 - a2 + new_a = self.set_bond_distance(new_a, atom.element) self.add_proton(atom, c+new_a) @@ -307,20 +310,20 @@ class Protonate: def tetrahedral(self, atom): - self.display('TETRAHEDRAL - %d bonded atoms'%(len(atom.bonded_atoms))) + self.display('TETRAHEDRAL - %d bonded atoms'%(len(atom.bonded_atoms))) rot_angle = math.radians(109.5) # sanity check # if atom.number_of_protons_to_add + len(atom.bonded_atoms) != 4: - # self.display 'Error: Attempting tetrahedral structure with %d bonds'%(atom.number_of_protons_to_add + + # self.display 'Error: Attempting tetrahedral structure with %d bonds'%(atom.number_of_protons_to_add + # len(atom.bonded_atoms)) - + c = vector(atom1 = atom) # 0 bonds if len(atom.bonded_atoms) == 0: pass - + # 1 bond if len(atom.bonded_atoms) == 1 and atom.number_of_protons_to_add > 0: # Add another atom with the right angle to the first one @@ -332,7 +335,7 @@ class Protonate: # 2 bonds if len(atom.bonded_atoms) == 2 and atom.number_of_protons_to_add > 0: - # Add another atom with the right angle to the first two + # Add another atom with the right angle to the first two a1 = vector(atom1 = atom, atom2 = atom.bonded_atoms[0]).rescale(1.0) a2 = vector(atom1 = atom, atom2 = atom.bonded_atoms[1]).rescale(1.0) @@ -351,16 +354,16 @@ class Protonate: new_a = -a1-a2-a3 new_a = self.set_bond_distance(new_a, atom.element) self.add_proton(atom, c+new_a) - + return def add_proton(self, atom, position): # Create the new proton - new_H = Source.atom.Atom() - new_H.setProperty(numb = None, - name = 'H%s'%atom.name[1:], - resName = atom.resName, + new_H = propka.atom.Atom() + new_H.setProperty(numb = None, + name = 'H%s'%atom.name[1:], + resName = atom.resName, chainID = atom.chainID, resNumb = atom.resNumb, x = round(position.x,3), # round of to three digimal points @@ -381,19 +384,19 @@ class Protonate: atom.bonded_atoms.append(new_H) atom.number_of_protons_to_add -=1 - atom.conformation_container.add_atom(new_H) + atom.conformation_container.add_atom(new_H) # update names of all protons on this atom - new_H.residue_label = "%-3s%4d%2s" % (new_H.name,new_H.resNumb, new_H.chainID) + new_H.residue_label = "%-3s%4d%2s" % (new_H.name,new_H.resNumb, new_H.chainID) no_protons = atom.count_bonded_elements('H') if no_protons > 1: i = 1 for proton in atom.get_bonded_elements('H'): proton.name = 'H%s%d'%(atom.name[1:],i) - proton.residue_label = "%-3s%4d%2s" % (proton.name,proton.resNumb, proton.chainID) + proton.residue_label = "%-3s%4d%2s" % (proton.name,proton.resNumb, proton.chainID) i+=1 - + self.display('added',new_H, 'to',atom) return @@ -407,7 +410,7 @@ class Protonate: a = a.rescale(d) return a - + def display(self,*text): if self.verbose: s = '' @@ -429,11 +432,11 @@ if __name__ == '__main__': print('Error: Could not find \"%s\"'%filename) sys.exit(1) - + p = Protonate() pdblist = pdb.readPDB(filename) my_protein = protein.Protein(pdblist,'test.pdb') - + p.remove_all_hydrogen_atoms_from_protein(my_protein) my_protein.writePDB('before_protonation.pdb') diff --git a/propka/run.py b/propka/run.py new file mode 100644 index 0000000..867a000 --- /dev/null +++ b/propka/run.py @@ -0,0 +1,15 @@ +# entry point for propka script + +import propka.lib, propka.molecular_container + +def main(): + """ + Reads in structure files, calculates pKa values, and prints pKa files + """ + # loading options, flaggs and arguments + options, pdbfiles = propka.lib.loadOptions() + + for pdbfile in pdbfiles: + my_molecule = propka.molecular_container.Molecular_container(pdbfile, options) + my_molecule.calculate_pka() + my_molecule.write_pka() diff --git a/Source/vector_algebra.py b/propka/vector_algebra.py similarity index 99% rename from Source/vector_algebra.py rename to propka/vector_algebra.py index 0603ce9..0b9f0bb 100644 --- a/Source/vector_algebra.py +++ b/propka/vector_algebra.py @@ -1,3 +1,5 @@ +from __future__ import division +from __future__ import print_function import math class vector: diff --git a/Source/version.py b/propka/version.py similarity index 80% rename from Source/version.py rename to propka/version.py index 5a5f603..2ec9410 100644 --- a/Source/version.py +++ b/propka/version.py @@ -1,8 +1,11 @@ +from __future__ import division +from __future__ import print_function import math -import Source.lib as lib import sys, os -import Source.calculations as calculations -import Source.parameters + +import propka.lib as lib +import propka.calculations as calculations +import propka.parameters class version: @@ -21,7 +24,7 @@ class version: def hydrogen_bond_interaction(self, group1, group2): return self.hydrogen_bond_interaction_model(group1, group2, self) - def calculateSideChainEnergy(self, distance, dpka_max, cutoff, weight, f_angle): + def calculateSideChainEnergy(self, distance, dpka_max, cutoff, weight, f_angle): return self.sidechain_interaction_model(distance, dpka_max, cutoff, f_angle) # weight is ignored in 3.0 Sep07 # coulomb @@ -47,38 +50,38 @@ class version: def setup_bonding(self, molecular_container): return self.prepare_bonds(self.parameters, molecular_container) - + class version_A(version): def __init__(self, parameters): - # set the calculation rutines used in this version + # set the calculation rutines used in this version version.__init__(self, parameters) # atom naming, bonding, and protonation - self.molecular_preparation_method = Source.calculations.setup_bonding_and_protonation - self.prepare_bonds = Source.calculations.setup_bonding + self.molecular_preparation_method = propka.calculations.setup_bonding_and_protonation + self.prepare_bonds = propka.calculations.setup_bonding # desolvation related methods self.desolvation_model = calculations.radial_volume_desolvation self.weight_pair_method = calculations.calculatePairWeight - + # side chain methods - self.sidechain_interaction_model = Source.calculations.HydrogenBondEnergy - self.hydrogen_bond_interaction_model = Source.calculations.hydrogen_bond_interaction + self.sidechain_interaction_model = propka.calculations.HydrogenBondEnergy + self.hydrogen_bond_interaction_model = propka.calculations.hydrogen_bond_interaction # colomb methods - self.electrostatic_interaction_model = Source.calculations.electrostatic_interaction - self.check_coulomb_pair_method = Source.calculations.checkCoulombPair - self.coulomb_interaction_model = Source.calculations.CoulombEnergy + self.electrostatic_interaction_model = propka.calculations.electrostatic_interaction + self.check_coulomb_pair_method = propka.calculations.checkCoulombPair + self.coulomb_interaction_model = propka.calculations.CoulombEnergy #backbone - self.backbone_interaction_model = Source.calculations.HydrogenBondEnergy - self.backbone_reorganisation_method = Source.calculations.BackBoneReorganization + self.backbone_interaction_model = propka.calculations.HydrogenBondEnergy + self.backbone_reorganisation_method = propka.calculations.BackBoneReorganization # exception methods - self.exception_check_method = Source.calculations.checkExceptions + self.exception_check_method = propka.calculations.checkExceptions return def get_hydrogen_bond_parameters(self, atom1, atom2): @@ -91,7 +94,7 @@ class version_A(version): if atom.group_type in self.parameters.backbone_CO_hydrogen_bond.keys(): [v,c1,c2] = self.parameters.backbone_CO_hydrogen_bond[atom.group_type] return [v,[c1,c2]] - + if backbone_atom.group_type == 'BBN': if atom.group_type in self.parameters.backbone_NH_hydrogen_bond.keys(): [v,c1,c2] = self.parameters.backbone_NH_hydrogen_bond[atom.group_type] @@ -104,7 +107,7 @@ class version_A(version): class simple_hb(version_A): def __init__(self, parameters): - # set the calculation rutines used in this version + # set the calculation rutines used in this version version_A.__init__(self, parameters) print('Using simple hb model') return @@ -116,23 +119,23 @@ class simple_hb(version_A): def get_backbone_hydrogen_bond_parameters(self, backbone_atom, atom): return self.parameters.hydrogen_bonds.get_value(backbone_atom.element, atom.element) - + class element_based_ligand_interactions(version_A): def __init__(self, parameters): - # set the calculation rutines used in this version + # set the calculation rutines used in this version version_A.__init__(self, parameters) print('Using detailed SC model!') return def get_hydrogen_bond_parameters(self, atom1, atom2): - if not 'hetatm' in [atom1.type, atom2.type]: + if not 'hetatm' in [atom1.type, atom2.type]: # this is a protein-protein interaction dpka_max = self.parameters.sidechain_interaction.get_value(atom1.group_type, atom2.group_type) cutoff = self.parameters.sidechain_cutoffs.get_value(atom1.group_type, atom2.group_type) return [dpka_max, cutoff] - + # at least one ligand atom is involved in this interaction # make sure that we are using the heavy atoms for finding paramters elements = [] @@ -140,17 +143,17 @@ class element_based_ligand_interactions(version_A): if a.element == 'H': elements.append(a.bonded_atoms[0].element) else: elements.append(a.element) - return self.parameters.hydrogen_bonds.get_value(elements[0], elements[1]) + return self.parameters.hydrogen_bonds.get_value(elements[0], elements[1]) + - def get_backbone_hydrogen_bond_parameters(self, backbone_atom, atom): - if atom.type == 'atom': + if atom.type == 'atom': # this is a backbone-protein interaction if backbone_atom.group_type == 'BBC' and\ atom.group_type in self.parameters.backbone_CO_hydrogen_bond.keys(): [v,c1,c2] = self.parameters.backbone_CO_hydrogen_bond[atom.group_type] return [v,[c1,c2]] - + if backbone_atom.group_type == 'BBN' and\ atom.group_type in self.parameters.backbone_NH_hydrogen_bond.keys(): [v,c1,c2] = self.parameters.backbone_NH_hydrogen_bond[atom.group_type] @@ -168,7 +171,7 @@ class element_based_ligand_interactions(version_A): print('Could not determine backbone interaction parameters for:', backbone_atom,atom) - return + return return None @@ -176,28 +179,28 @@ class element_based_ligand_interactions(version_A): class propka30(version): def __init__(self, parameters): - # set the calculation rutines used in this version + # set the calculation rutines used in this version version.__init__(self, parameters) # atom naming, bonding, and protonation - self.molecular_preparation_method = Source.calculations.setup_bonding_and_protonation_30_style + self.molecular_preparation_method = propka.calculations.setup_bonding_and_protonation_30_style # desolvation related methods self.desolvation_model = calculations.radial_volume_desolvation self.weight_pair_method = calculations.calculatePairWeight - + # side chain methods - self.sidechain_interaction_model = Source.calculations.HydrogenBondEnergy + self.sidechain_interaction_model = propka.calculations.HydrogenBondEnergy # colomb methods - self.check_coulomb_pair_method = Source.calculations.checkCoulombPair - self.coulomb_interaction_model = Source.calculations.CoulombEnergy + self.check_coulomb_pair_method = propka.calculations.checkCoulombPair + self.coulomb_interaction_model = propka.calculations.CoulombEnergy #backbone - self.backbone_reorganisation_method = Source.calculations.BackBoneReorganization + self.backbone_reorganisation_method = propka.calculations.BackBoneReorganization # exception methods - self.exception_check_method = Source.calculations.checkExceptions + self.exception_check_method = propka.calculations.checkExceptions return diff --git a/scripts/propka31.py b/scripts/propka31.py new file mode 100755 index 0000000..a35426b --- /dev/null +++ b/scripts/propka31.py @@ -0,0 +1,31 @@ +#!/usr/bin/env python + +# This is the original propka script. However, this distribute-based +# installation moved the main() function into propka.run.main and just +# generates a script called propka31 from the setup.py installation +# script. You should not need to use this script. +# +# (Also note that there can be import problems because the script name +# is the same as the module name; that's why the new script is called +# propka31.) + +import propka.lib, propka.molecular_container + +def main(): + """ + Reads in structure files, calculates pKa values, and prints pKa files + """ + # loading options, flaggs and arguments + options, pdbfiles = propka.lib.loadOptions() + + for pdbfile in pdbfiles: + my_molecule = propka.molecular_container.Molecular_container(pdbfile, options) + my_molecule.calculate_pka() + my_molecule.write_pka() + + +if __name__ == '__main__': + #import cProfile + #cProfile.run('main()',sort=1) + main() + diff --git a/setup.py b/setup.py new file mode 100644 index 0000000..1bee359 --- /dev/null +++ b/setup.py @@ -0,0 +1,51 @@ +# PROPKA 3.1 +# +# +# setuptools installation of PROPKA 3.1 + +import ez_setup +ez_setup.use_setuptools() +from setuptools import setup, find_packages + +VERSION = "3.1" + +setup(name="PROPKA", + version=VERSION, + description="Heuristic pKa calculations with ligands", + long_description=""" +PROPKA predicts the pKa values of ionizable groups in proteins (version 3.0) and +protein-ligand complexes (version 3.1) based on the 3D structure. + +For proteins without ligands both version should produce the same result. + +The method is described in the following papers, which you should cite +in publications: + +* Sondergaard, Chresten R., Mats HM Olsson, Michal Rostkowski, and Jan + H. Jensen. "Improved Treatment of Ligands and Coupling Effects in + Empirical Calculation and Rationalization of pKa Values." Journal of + Chemical Theory and Computation 7, no. 7 (2011): 2284-2295. + +* Olsson, Mats HM, Chresten R. Sondergaard, Michal Rostkowski, and Jan + H. Jensen. "PROPKA3: consistent treatment of internal and surface + residues in empirical pKa predictions." Journal of Chemical Theory + and Computation 7, no. 2 (2011): 525-537. + +See http://propka.ki.ku.dk/ for the PROPKA web server, +using the tutorial http://propka.ki.ku.dk/~luca/wiki/index.php/PROPKA_3.1_Tutorial . +""", + author="Jan H. Jensen", + author_email="jhjensen@chem.ku.dk", + license="", + url="http://propka.ki.ku.dk/", + keywords="science", + packages=find_packages(exclude=['scripts']), + package_data = {'propka': ['*.dat', '*.cfg']}, + #scripts = ["scripts/propka31.py"], # use entry point below + entry_points = { + 'console_scripts': [ + 'propka31 = propka.run:main', + ], + }, + zip_safe=True, +)