diff --git a/.pymol/startup/Draw_Protein_Dimensions.py b/.pymol/startup/Draw_Protein_Dimensions.py new file mode 100644 index 0000000..3014a49 --- /dev/null +++ b/.pymol/startup/Draw_Protein_Dimensions.py @@ -0,0 +1,362 @@ + +''' +Calculate and display the dimensions of a protein. + + +This is a first version, please use at your own risk! + +REQUIREMENTS + +numpy (http://numpy.scipy.org) that should be built into the newers versions of Pymol + +(c) Pablo Guardado Calvo + +Based on "inertia_tensor.py" (c) 2010 by Mateusz Maciejewski + +License: MIT + +''' + +from __future__ import print_function + +__author__ = 'Pablo Guardado Calvo' +__version__ = '0.1' +__email__ = 'pablo.guardado (at) gmail.com' +__date__ = '13/08/2015' + +########################################################################################################################################################### +# USAGE +# +# The idea behing this script is to calculate an aproximate minimal bounding box to extract the cell dimensions of a protein. To calculate the minimal bounding +# is not trivial and usually the Axis Aligned Bounding Box (AABB) does not show up the real dimensions of the protein. This script calculates the inertia tensor +# of the object, extract the eigenvalues and use them to rotate the molecule (using as rotation matrix the transpose of the eigenvalues matrix). The result is that +# the molecule is oriented with the inertia axis aligned with the cartesian axis. A new Bounding Box is calculated that is called Inertia Axis Aligned Bounding Box +#(IABB), whose volume is always lower than AABB volume, and in many cases will correspond with the lowest volume. Of course, maybe it exists another Bounding Box +# with a lower volume (the minimal Bounding Box). +# +# As always with these type of things, you have to use at your own risk. I did not try all the possible combinations, but if you find a bug, do +# not hesitate to contact me (pablo.guardado (at) gmail.com) or try to modify the code for yourself to correct it. +# +# To load the script just type: +# +# run path-to-the-script/Draw_Protein_Dimensions.py +# +# or if you want something more permanent add the previous line to your .pymolrc file +# +# The script works just typing: +# +# draw_Protein_Dimensions selection +# +# This will draw the cell dimensions of your selection based on a IABB. It also generates the IABB box and the inertia axis, you just need to do "show cgo" to display them. +# +# You could also try: +# +# draw_BB selection +# +# This will draw the AABB and IABB boxes with their cell dimensions and show in the command line their volumes, you can compare both of them. +############################################################################################################################################################ + +from pymol import cmd, cgo +from pymol.cgo import * +import numpy +from random import randint + +def matriz_inercia(selection): + ''' + DESCRIPTION + + The method calculates the mass center, the inertia tensor and the eigenvalues and eigenvectors + for a given selection. Mostly taken from inertia_tensor.py + ''' + + model = cmd.get_model(selection) + totmass = 0.0 + x,y,z = 0,0,0 + for a in model.atom: + m = a.get_mass() + x += a.coord[0]*m + y += a.coord[1]*m + z += a.coord[2]*m + totmass += m + global cM + cM = numpy.array([x/totmass, y/totmass, z/totmass]) + + + I = [] + for index in range(9): + I.append(0) + + for a in model.atom: + temp_x, temp_y, temp_z = a.coord[0], a.coord[1], a.coord[2] + temp_x -= x + temp_y -= y + temp_z -= z + + I[0] += a.get_mass() * (temp_y**2 + temp_z**2) + I[1] -= a.get_mass() * temp_x * temp_y + I[2] -= a.get_mass() * temp_x * temp_z + I[3] -= a.get_mass() * temp_x * temp_y + I[4] += a.get_mass() * (temp_x**2 + temp_z**2) + I[5] -= a.get_mass() * temp_y * temp_z + I[6] -= a.get_mass() * temp_x * temp_z + I[7] -= a.get_mass() * temp_y * temp_z + I[8] += a.get_mass() * (temp_x**2 + temp_y**2) + + global tensor + tensor = numpy.array([(I[0:3]), (I[3:6]), (I[6:9])]) + + global autoval, autovect, ord_autoval, ord_autovect + autoval, autovect = numpy.linalg.eig(tensor) + auto_ord = numpy.argsort(autoval) + ord_autoval = autoval[auto_ord] + ord_autovect_complete = autovect[:, auto_ord].T + ord_autovect = numpy.around(ord_autovect_complete, 3) + + return ord_autoval + + +def draw_inertia_axis(selection): + ''' + DESCRIPTION + + This method draw the inertia axis calculated with the method matriz_inercia. + ''' + + matriz_inercia(selection) + axis1 = ord_autovect[0] + x1, y1, z1 = cM[0], cM[1], cM[2] + x2, y2, z2 = cM[0]+50*axis1[0], cM[1]+50*axis1[1], cM[2]+50*axis1[2] + eje1 = [cgo.CYLINDER, x1, y1, z1, x2, y2, z2, 0.6, 1, 0, 0, 1, 0, 0, 0.0] + cmd.load_cgo(eje1, 'Inertia_Axis1') + axis2 = ord_autovect[1] + x3, y3, z3 = cM[0]+40*axis2[0], cM[1]+40*axis2[1], cM[2]+40*axis2[2] + eje1 = [cgo.CYLINDER, x1, y1, z1, x3, y3, z3, 0.6, 1, 0.5, 0, 1, 0.5, 0, 0.0] + cmd.load_cgo(eje1, 'Inertia_Axis2') + axis4 = ord_autovect[2] + x4, y4, z4 = cM[0]+30*axis4[0], cM[1]+30*axis4[1], cM[2]+30*axis4[2] + eje1 = [cgo.CYLINDER, x1, y1, z1, x4, y4, z4, 0.6, 1, 1, 0, 1, 1, 0, 0.0] + cmd.load_cgo(eje1, 'Inertia_Axis3') + +def translacion_cM(selection): + ''' + DESCRIPTION + + Translate the center of mass of the molecule to the origin. + ''' + model = cmd.get_model(selection) + totmass = 0.0 + x,y,z = 0,0,0 + for a in model.atom: + m = a.get_mass() + x += a.coord[0]*m + y += a.coord[1]*m + z += a.coord[2]*m + totmass += m + cM = numpy.array([x/totmass, y/totmass, z/totmass]) + trans_array = ([1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, -cM[0], -cM[1], -cM[2], 1]) + model_trans = cmd.transform_selection(selection, trans_array) + + +def rotacion_orig(selection): + ''' + DESCRIPTION + + Find the proper rotation matrix, i.e. the transpose of the matrix formed by the eigenvectors of the inertia tensor + ''' + + translacion_cM(selection) + matriz_inercia(selection) + global transf, transf_array, ord_autovect_array, transf_array_print + ord_autovect_array = numpy.array([[ord_autovect[0][0], ord_autovect[0][1], ord_autovect[0][2]], + [ord_autovect[1][0], ord_autovect[1][1], ord_autovect[1][2]], + [ord_autovect[2][0], ord_autovect[2][1], ord_autovect[2][2]]]) + if numpy.linalg.det(ord_autovect_array) == -1: + ord_autovect_array = numpy.array([[ord_autovect[2][0], ord_autovect[2][1], ord_autovect[2][2]], + [ord_autovect[1][0], ord_autovect[1][1], ord_autovect[1][2]], + [ord_autovect[0][0], ord_autovect[0][1], ord_autovect[0][2]]]) + transf = numpy.transpose(ord_autovect_array) + transf_array = numpy.array([transf[0][0], transf[0][1], transf[0][2], 0, + transf[1][0], transf[1][1], transf[1][2], 0, + transf[2][0], transf[2][1], transf[2][2], 0, + 0, 0, 0, 1]) + + +def transformar(selection): + ''' + DESCRIPTION + + Rotate the molecule and draw the inertia axis. + ''' + + rotacion_orig(selection) + model_rot = cmd.transform_selection(selection, transf_array, homogenous=0, transpose=1); + draw_inertia_axis(selection) + + + +def draw_AABB(selection): + """ + DESCRIPTION + For a given selection, draw the Axes Aligned bounding box around it without padding. Code taken and modified from DrawBoundingBox.py. + + """ + + AA_original = selection + "_original" + model_orig = cmd.create(AA_original, selection) + + ([min_X, min_Y, min_Z],[max_X, max_Y, max_Z]) = cmd.get_extent(AA_original) + + print("The Axis Aligned Bounding Box (AABB) dimensions are (%.2f, %.2f, %.2f)" % (max_X-min_X, max_Y-min_Y, max_Z-min_Z)) + print("The Axis Aligned Bounding Box (AABB) volume is %.2f A3" % ((max_X-min_X)*(max_Y-min_Y)*(max_Z-min_Z))) + + + min_X = min_X + min_Y = min_Y + min_Z = min_Z + max_X = max_X + max_Y = max_Y + max_Z = max_Z + + boundingBox = [ + LINEWIDTH, float(2), + + BEGIN, LINES, + COLOR, float(1), float(1), float(0), + + VERTEX, min_X, min_Y, min_Z, + VERTEX, min_X, min_Y, max_Z, + VERTEX, min_X, max_Y, min_Z, + VERTEX, min_X, max_Y, max_Z, + VERTEX, max_X, min_Y, min_Z, + VERTEX, max_X, min_Y, max_Z, + VERTEX, max_X, max_Y, min_Z, + VERTEX, max_X, max_Y, max_Z, + VERTEX, min_X, min_Y, min_Z, + VERTEX, max_X, min_Y, min_Z, + VERTEX, min_X, max_Y, min_Z, + VERTEX, max_X, max_Y, min_Z, + VERTEX, min_X, max_Y, max_Z, + VERTEX, max_X, max_Y, max_Z, + VERTEX, min_X, min_Y, max_Z, + VERTEX, max_X, min_Y, max_Z, + VERTEX, min_X, min_Y, min_Z, + VERTEX, min_X, max_Y, min_Z, + VERTEX, max_X, min_Y, min_Z, + VERTEX, max_X, max_Y, min_Z, + VERTEX, min_X, min_Y, max_Z, + VERTEX, min_X, max_Y, max_Z, + VERTEX, max_X, min_Y, max_Z, + VERTEX, max_X, max_Y, max_Z, + + END + ] + + p0 = '_0' + str(randint(0, 100)) + p1 = '_1' + str(randint(0, 100)) + p2 = '_2' + str(randint(0, 100)) + p3 = '_3' + str(randint(0, 100)) + cmd.pseudoatom (pos=[min_X, min_Y, min_Z], object=p0) + cmd.pseudoatom (pos=[min_X, min_Y, max_Z], object=p1) + cmd.pseudoatom (pos=[min_X, max_Y, min_Z], object=p2) + cmd.pseudoatom (pos=[max_X, min_Y, min_Z], object=p3) + cmd.distance(None, p0, p3) + cmd.distance(None, p0, p2) + cmd.distance(None, p0, p1) + cmd.hide("nonbonded") + + boxName = "box_AABB_" + str(randint(0, 100)) + cmd.load_cgo(boundingBox,boxName) + return boxName + + +def draw_IABB(selection): + """ + DESCRIPTION + For a given selection, draw the Inertia Axes Aligned bounding box around it without padding. Code taken and modified from DrawBoundingBox.py. + + """ + + transformar(selection) + + ([minX, minY, minZ],[maxX, maxY, maxZ]) = cmd.get_extent(selection) + + print("The Inertia Axis Aligned Bounding Box (IABB) dimensions are (%.2f, %.2f, %.2f)" % (maxX-minX, maxY-minY, maxZ-minZ)) + print("The Inertia Axis Aligned Bounding Box (IABB) volume is %.2f A3" % ((maxX-minX)*(maxY-minY)*(maxZ-minZ))) + + + minX = minX + minY = minY + minZ = minZ + maxX = maxX + maxY = maxY + maxZ = maxZ + + boundingBox = [ + LINEWIDTH, float(2), + + BEGIN, LINES, + COLOR, float(1), float(0), float(0), + + VERTEX, minX, minY, minZ, + VERTEX, minX, minY, maxZ, + VERTEX, minX, maxY, minZ, + VERTEX, minX, maxY, maxZ, + VERTEX, maxX, minY, minZ, + VERTEX, maxX, minY, maxZ, + VERTEX, maxX, maxY, minZ, + VERTEX, maxX, maxY, maxZ, + VERTEX, minX, minY, minZ, + VERTEX, maxX, minY, minZ, + VERTEX, minX, maxY, minZ, + VERTEX, maxX, maxY, minZ, + VERTEX, minX, maxY, maxZ, + VERTEX, maxX, maxY, maxZ, + VERTEX, minX, minY, maxZ, + VERTEX, maxX, minY, maxZ, + VERTEX, minX, minY, minZ, + VERTEX, minX, maxY, minZ, + VERTEX, maxX, minY, minZ, + VERTEX, maxX, maxY, minZ, + VERTEX, minX, minY, maxZ, + VERTEX, minX, maxY, maxZ, + VERTEX, maxX, minY, maxZ, + VERTEX, maxX, maxY, maxZ, + + END + ] + + p4 = '_4' + str(randint(0, 100)) + p5 = '_5' + str(randint(0, 100)) + p6 = '_6' + str(randint(0, 100)) + p7 = '_7' + str(randint(0, 100)) + cmd.pseudoatom (pos=[minX, minY, minZ], object=p4) + cmd.pseudoatom (pos=[minX, minY, maxZ], object=p5) + cmd.pseudoatom (pos=[minX, maxY, minZ], object=p6) + cmd.pseudoatom (pos=[maxX, minY, minZ], object=p7) + cmd.distance(None, p4, p7) + cmd.distance(None, p4, p6) + cmd.distance(None, p4, p5) + cmd.hide("nonbonded") + + boxName = "box_IABB_" + str(randint(0, 100)) + cmd.load_cgo(boundingBox,boxName) + return boxName + + +def draw_BB(selection): + draw_AABB(selection) + draw_IABB(selection) + +def draw_Protein_Dimensions(selection): + draw_IABB(selection) + cmd.hide("cgo") + +cmd.extend ("draw_Protein_Dimensions", draw_Protein_Dimensions) +cmd.extend ("draw_BB", draw_BB) + + + + + + + diff --git a/.pymol/startup/Draw_Protein_Dimensions.pyc b/.pymol/startup/Draw_Protein_Dimensions.pyc new file mode 100644 index 0000000..5479534 Binary files /dev/null and b/.pymol/startup/Draw_Protein_Dimensions.pyc differ diff --git a/.pymol/startup/axes.py b/.pymol/startup/axes.py new file mode 100644 index 0000000..bc2e0a4 --- /dev/null +++ b/.pymol/startup/axes.py @@ -0,0 +1,68 @@ +from pymol import cmd +from chempy import cpv + +class PutCenterCallback(object): + prev_v = None + + def __init__(self, name, corner=0): + self.name = name + self.corner = corner + self.cb_name = cmd.get_unused_name('_cb') + + def load(self): + cmd.load_callback(self, self.cb_name) + + def __call__(self): + if self.name not in cmd.get_names('objects'): + import threading + threading.Thread(None, cmd.delete, args=(self.cb_name,)).start() + return + + v = cmd.get_view() + if v == self.prev_v: + return + self.prev_v = v + + t = v[12:15] + + if self.corner: + vp = cmd.get_viewport() + R_mc = [v[0:3], v[3:6], v[6:9]] + off_c = [0.15 * v[11] * vp[0] / vp[1], 0.15 * v[11], 0.0] + if self.corner in [2,3]: + off_c[0] *= -1 + if self.corner in [3,4]: + off_c[1] *= -1 + off_m = cpv.transform(R_mc, off_c) + t = cpv.add(t, off_m) + + z = -v[11] / 30.0 + m = [z, 0, 0, t[0] / z, 0, z, 0, t[1] / z, 0, 0, z, t[2] / z, 0, 0, 0, 1] + cmd.set_object_ttt(self.name, m, homogenous=1) + +def axes(name='axes'): + ''' +DESCRIPTION + + Puts coordinate axes to the lower left corner of the viewport. + ''' + from pymol import cgo + + cmd.set('auto_zoom', 0) + + w = 0.06 # cylinder width + l = 0.75 # cylinder length + h = 0.25 # cone hight + d = w * 1.618 # cone base diameter + + obj = [cgo.CYLINDER, 0.0, 0.0, 0.0, l, 0.0, 0.0, w, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, + cgo.CYLINDER, 0.0, 0.0, 0.0, 0.0, l, 0.0, w, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, + cgo.CYLINDER, 0.0, 0.0, 0.0, 0.0, 0.0, l, w, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, + cgo.CONE, l, 0.0, 0.0, h+l, 0.0, 0.0, d, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, + cgo.CONE, 0.0, l, 0.0, 0.0, h+l, 0.0, d, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, + cgo.CONE, 0.0, 0.0, l, 0.0, 0.0, h+l, d, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 1.0] + + PutCenterCallback(name, 1).load() + cmd.load_cgo(obj, name) + +cmd.extend('axes', axes) diff --git a/.pymol/startup/axes.pyc b/.pymol/startup/axes.pyc new file mode 100644 index 0000000..79df320 Binary files /dev/null and b/.pymol/startup/axes.pyc differ diff --git a/.pymol/startup/center_of_mass.py b/.pymol/startup/center_of_mass.py new file mode 100644 index 0000000..e9b103f --- /dev/null +++ b/.pymol/startup/center_of_mass.py @@ -0,0 +1,86 @@ +''' +See more here: http://www.pymolwiki.org/index.php/center_of_mass + +DESCRIPTION + + Places a pseudoatom at the center of mass + + Author: Sean Law + Michigan State University + slaw (at) msu . edu + +SEE ALSO + + pseudoatom, get_com +''' + +from __future__ import print_function +from pymol import cmd + + +def com(selection, state=None, mass=None, object=None, quiet=1, **kwargs): + quiet = int(quiet) + if (object == None): + try: + object = cmd.get_legal_name(selection) + object = cmd.get_unused_name(object + "_COM", 0) + except AttributeError: + object = 'COM' + cmd.delete(object) + + if (state != None): + x, y, z = get_com(selection, mass=mass, quiet=quiet) + if not quiet: + print("%f %f %f" % (x, y, z)) + cmd.pseudoatom(object, pos=[x, y, z], **kwargs) + cmd.show("spheres", object) + else: + for i in range(cmd.count_states()): + x, y, z = get_com(selection, mass=mass, state=i + 1, quiet=quiet) + if not quiet: + print("State %d:%f %f %f" % (i + 1, x, y, z)) + cmd.pseudoatom(object, pos=[x, y, z], state=i + 1, **kwargs) + cmd.show("spheres", 'last ' + object) + +cmd.extend("com", com) + + +def get_com(selection, state=1, mass=None, quiet=1): + """ + DESCRIPTION + + Calculates the center of mass + + Author: Sean Law + Michigan State University + slaw (at) msu . edu + """ + quiet = int(quiet) + + totmass = 0.0 + if mass != None and not quiet: + print("Calculating mass-weighted COM") + + state = int(state) + model = cmd.get_model(selection, state) + x, y, z = 0, 0, 0 + for a in model.atom: + if (mass != None): + m = a.get_mass() + x += a.coord[0] * m + y += a.coord[1] * m + z += a.coord[2] * m + totmass += m + else: + x += a.coord[0] + y += a.coord[1] + z += a.coord[2] + + if (mass != None): + return x / totmass, y / totmass, z / totmass + else: + return x / len(model.atom), y / len(model.atom), z / len(model.atom) + +cmd.extend("get_com", get_com) + +# vi:expandtab:sw=3 diff --git a/.pymol/startup/center_of_mass.pyc b/.pymol/startup/center_of_mass.pyc new file mode 100644 index 0000000..a853fa8 Binary files /dev/null and b/.pymol/startup/center_of_mass.pyc differ diff --git a/.pymol/startup/color_h.py b/.pymol/startup/color_h.py new file mode 100644 index 0000000..160e866 --- /dev/null +++ b/.pymol/startup/color_h.py @@ -0,0 +1,128 @@ +# color_h +# ------- + +# PyMOL command to color protein molecules according to the Eisenberg hydrophobicity scale + +# +# Source: http://us.expasy.org/tools/pscale/Hphob.Eisenberg.html +# Amino acid scale: Normalized consensus hydrophobicity scale +# Author(s): Eisenberg D., Schwarz E., Komarony M., Wall R. +# Reference: J. Mol. Biol. 179:125-142 (1984) +# +# Amino acid scale values: +# +# Ala: 0.620 +# Arg: -2.530 +# Asn: -0.780 +# Asp: -0.900 +# Cys: 0.290 +# Gln: -0.850 +# Glu: -0.740 +# Gly: 0.480 +# His: -0.400 +# Ile: 1.380 +# Leu: 1.060 +# Lys: -1.500 +# Met: 0.640 +# Phe: 1.190 +# Pro: 0.120 +# Ser: -0.180 +# Thr: -0.050 +# Trp: 0.810 +# Tyr: 0.260 +# Val: 1.080 +# +# Usage: +# color_h (selection) +# +from pymol import cmd + +def color_h(selection='all'): + s = str(selection) + print s + cmd.set_color('color_ile',[0.996,0.062,0.062]) + cmd.set_color('color_phe',[0.996,0.109,0.109]) + cmd.set_color('color_val',[0.992,0.156,0.156]) + cmd.set_color('color_leu',[0.992,0.207,0.207]) + cmd.set_color('color_trp',[0.992,0.254,0.254]) + cmd.set_color('color_met',[0.988,0.301,0.301]) + cmd.set_color('color_ala',[0.988,0.348,0.348]) + cmd.set_color('color_gly',[0.984,0.394,0.394]) + cmd.set_color('color_cys',[0.984,0.445,0.445]) + cmd.set_color('color_tyr',[0.984,0.492,0.492]) + cmd.set_color('color_pro',[0.980,0.539,0.539]) + cmd.set_color('color_thr',[0.980,0.586,0.586]) + cmd.set_color('color_ser',[0.980,0.637,0.637]) + cmd.set_color('color_his',[0.977,0.684,0.684]) + cmd.set_color('color_glu',[0.977,0.730,0.730]) + cmd.set_color('color_asn',[0.973,0.777,0.777]) + cmd.set_color('color_gln',[0.973,0.824,0.824]) + cmd.set_color('color_asp',[0.973,0.875,0.875]) + cmd.set_color('color_lys',[0.899,0.922,0.922]) + cmd.set_color('color_arg',[0.899,0.969,0.969]) + cmd.color("color_ile","("+s+" and resn ile)") + cmd.color("color_phe","("+s+" and resn phe)") + cmd.color("color_val","("+s+" and resn val)") + cmd.color("color_leu","("+s+" and resn leu)") + cmd.color("color_trp","("+s+" and resn trp)") + cmd.color("color_met","("+s+" and resn met)") + cmd.color("color_ala","("+s+" and resn ala)") + cmd.color("color_gly","("+s+" and resn gly)") + cmd.color("color_cys","("+s+" and resn cys)") + cmd.color("color_tyr","("+s+" and resn tyr)") + cmd.color("color_pro","("+s+" and resn pro)") + cmd.color("color_thr","("+s+" and resn thr)") + cmd.color("color_ser","("+s+" and resn ser)") + cmd.color("color_his","("+s+" and resn his)") + cmd.color("color_glu","("+s+" and resn glu)") + cmd.color("color_asn","("+s+" and resn asn)") + cmd.color("color_gln","("+s+" and resn gln)") + cmd.color("color_asp","("+s+" and resn asp)") + cmd.color("color_lys","("+s+" and resn lys)") + cmd.color("color_arg","("+s+" and resn arg)") +cmd.extend('color_h',color_h) + +def color_h2(selection='all'): + s = str(selection) + print s + cmd.set_color("color_ile2",[0.938,1,0.938]) + cmd.set_color("color_phe2",[0.891,1,0.891]) + cmd.set_color("color_val2",[0.844,1,0.844]) + cmd.set_color("color_leu2",[0.793,1,0.793]) + cmd.set_color("color_trp2",[0.746,1,0.746]) + cmd.set_color("color_met2",[0.699,1,0.699]) + cmd.set_color("color_ala2",[0.652,1,0.652]) + cmd.set_color("color_gly2",[0.606,1,0.606]) + cmd.set_color("color_cys2",[0.555,1,0.555]) + cmd.set_color("color_tyr2",[0.508,1,0.508]) + cmd.set_color("color_pro2",[0.461,1,0.461]) + cmd.set_color("color_thr2",[0.414,1,0.414]) + cmd.set_color("color_ser2",[0.363,1,0.363]) + cmd.set_color("color_his2",[0.316,1,0.316]) + cmd.set_color("color_glu2",[0.27,1,0.27]) + cmd.set_color("color_asn2",[0.223,1,0.223]) + cmd.set_color("color_gln2",[0.176,1,0.176]) + cmd.set_color("color_asp2",[0.125,1,0.125]) + cmd.set_color("color_lys2",[0.078,1,0.078]) + cmd.set_color("color_arg2",[0.031,1,0.031]) + cmd.color("color_ile2","("+s+" and resn ile)") + cmd.color("color_phe2","("+s+" and resn phe)") + cmd.color("color_val2","("+s+" and resn val)") + cmd.color("color_leu2","("+s+" and resn leu)") + cmd.color("color_trp2","("+s+" and resn trp)") + cmd.color("color_met2","("+s+" and resn met)") + cmd.color("color_ala2","("+s+" and resn ala)") + cmd.color("color_gly2","("+s+" and resn gly)") + cmd.color("color_cys2","("+s+" and resn cys)") + cmd.color("color_tyr2","("+s+" and resn tyr)") + cmd.color("color_pro2","("+s+" and resn pro)") + cmd.color("color_thr2","("+s+" and resn thr)") + cmd.color("color_ser2","("+s+" and resn ser)") + cmd.color("color_his2","("+s+" and resn his)") + cmd.color("color_glu2","("+s+" and resn glu)") + cmd.color("color_asn2","("+s+" and resn asn)") + cmd.color("color_gln2","("+s+" and resn gln)") + cmd.color("color_asp2","("+s+" and resn asp)") + cmd.color("color_lys2","("+s+" and resn lys)") + cmd.color("color_arg2","("+s+" and resn arg)") +cmd.extend('color_h2',color_h2) \ No newline at end of file diff --git a/.pymol/startup/color_h.pyc b/.pymol/startup/color_h.pyc new file mode 100644 index 0000000..39f3290 Binary files /dev/null and b/.pymol/startup/color_h.pyc differ diff --git a/.pymol/startup/dssp_stride.py b/.pymol/startup/dssp_stride.py new file mode 100644 index 0000000..20c2def --- /dev/null +++ b/.pymol/startup/dssp_stride.py @@ -0,0 +1,981 @@ +""" 2011_04_11: Hongbo Zhu + PyMOL plugin for running DSSP and display SS using different colors. + The DSSP codes for secondary structure are: + + H - Alpha helix (4-12) + G - 3-10 helix + I - pi helix + E - Strand + B - Isolated beta-bridge residue + T - Turn + S - Bend + - - None + 2011_04_24: add stride + STRIDE code (taken from stride.doc) is nearly the same as DSSP + H Alpha helix + G 3-10 helix + I PI-helix + E Extended conformation + B or b Isolated bridge + T Turn + C Coil (none of the above) +""" + +# Copyright Notice +# ================ +# +# The PyMOL Plugin source code in this file is copyrighted, but you can +# freely use and copy it as long as you don't change or remove any of +# the copyright notices. +# +# ---------------------------------------------------------------------- +# This PyMOL Plugin is Copyright (C) 2011 by +# Hongbo Zhu +# +# All Rights Reserved +# +# Permission to use, copy, modify, distribute, and distribute modified +# versions of this software and its documentation for any purpose and +# without fee is hereby granted, provided that the above copyright +# notice appear in all copies and that both the copyright notice and +# this permission notice appear in supporting documentation, and that +# the name(s) of the author(s) not be used in advertising or publicity +# pertaining to distribution of the software without specific, written +# prior permission. +# +# THE AUTHOR(S) DISCLAIM ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, +# INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS. IN +# NO EVENT SHALL THE AUTHOR(S) BE LIABLE FOR ANY SPECIAL, INDIRECT OR +# CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF +# USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR +# OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR +# PERFORMANCE OF THIS SOFTWARE. +# ---------------------------------------------------------------------- + +# python lib +import os +import sys +import platform +if sys.version_info >= (2,4): + import subprocess # subprocess is introduced in python 2.4 +import math +import random +import tempfile +import tkSimpleDialog +import tkMessageBox +import tkFileDialog +import tkColorChooser + +import Tkinter + +# pymol lib +try: + from pymol import cmd + from pymol.cgo import * +except ImportError: + print 'Warning: pymol library cmd not found.' + sys.exit(1) + +# external lib +try: + import Pmw +except ImportError: + print 'Warning: failed to import Pmw. Exit ...' + sys.exit(1) + +VERBOSE = True + +################# +## here we go +################# +def __init__(self): + """ DSSP and Stride plugin for PyMol + """ + self.menuBar.addmenuitem('Plugin', 'command', + 'DSSP & Stride', label = 'DSSP & Stride', + command = lambda s=self : DSSPPlugin(s)) + + +################# +## GUI related +################# +class DSSPPlugin: + + def __init__(self, app): + + self.parent = app.root + self.dialog = Pmw.Dialog(self.parent, + buttons = ('Run DSSP', + 'Run Stride', + 'Update Color', + 'Update ss', + 'Exit'), + title = 'DSSP and Stride Plugin for PyMOL', + command = self.execute) + Pmw.setbusycursorattributes(self.dialog.component('hull')) + + # parameters used by DSSP + self.pymol_sel = Tkinter.StringVar() + self.dssp_bin = Tkinter.StringVar() + self.stride_bin = Tkinter.StringVar() + self.dssp_rlt_dict = {} + self.stride_rlt_dict = {} + self.ss_asgn_prog = None # which program is used to assign ss + + self.sel_obj_list = [] # there may be more than one seletion or object + # defined by self.pymol_sel + # treat each selection and object separately + if 'DSSP_BIN' not in os.environ and 'PYMOL_GIT_MOD' in os.environ: + if sys.platform.startswith('linux') and platform.machine() == 'x86_32': + initialdir_dssp = os.path.join(os.environ['PYMOL_GIT_MOD'],"DSSP","i86Linux2","dssp-2") + os.environ['DSSP_BIN'] = initialdir_dssp + elif sys.platform.startswith('linux') and platform.machine() == 'x86_64': + initialdir_dssp = os.path.join(os.environ['PYMOL_GIT_MOD'],"DSSP","ia64Linux2","dssp-2") + os.environ['DSSP_BIN'] = initialdir_dssp + elif sys.platform.startswith('win'): + initialdir_dssp = os.path.join(os.environ['PYMOL_GIT_MOD'],"DSSP","win32","dssp.exe") + os.environ['DSSP_BIN'] = initialdir_dssp + else: + pass + if 'DSSP_BIN' in os.environ: + if VERBOSE: print 'Found DSSP_BIN in environmental variables', os.environ['DSSP_BIN'] + self.dssp_bin.set(os.environ['DSSP_BIN']) + else: + if VERBOSE: print 'DSSP_BIN not found in environmental variables.' + self.dssp_bin.set('') + if 'STRIDE_BIN' not in os.environ and 'PYMOL_GIT_MOD' in os.environ: + if sys.platform.startswith('linux') and platform.machine() == 'x86_32': + initialdir_stride = os.path.join(os.environ['PYMOL_GIT_MOD'],"Stride","i86Linux2","stride") + os.environ['STRIDE_BIN'] = initialdir_stride + elif sys.platform.startswith('linux') and platform.machine() == 'x86_64': + initialdir_stride = os.path.join(os.environ['PYMOL_GIT_MOD'],"Stride","ia64Linux2","stride") + os.environ['STRIDE_BIN'] = initialdir_stride + elif sys.platform.startswith('win'): + initialdir_stride = os.path.join(os.environ['PYMOL_GIT_MOD'],"Stride","win32","stride.exe") + os.environ['STRIDE_BIN'] = initialdir_stride + else: + pass + if 'STRIDE_BIN' in os.environ: + if VERBOSE: print 'Found STRIDE_BIN in environmental variables', os.environ['STRIDE_BIN'] + self.stride_bin.set(os.environ['STRIDE_BIN']) + else: + if VERBOSE: print 'STRIDE_BIN not found in environmental variables.' + self.stride_bin.set('') + + # DSSP visualization color + # - H Alpha helix (4-12) + # - G 3-10 helix + # - I pi helix + # + # - E Extended strand + # - B Isolated beta-bridge residue + # + # - T Turn + # - S Bend + # - - None + # STRIDE does not have S and None, but it has two more + # codes: 'b' same as 'B' and 'C' for coil + self.DSSP_SSE_list = ['H','G','I', 'E','B','T','S','-'] # for the record of the order + self.STRIDE_SSE_list = ['H','G','I', 'E','B','b','T','C'] # for the record of the order + self.SSE_map = {'H':'H', + 'G':'H', + 'I':'H', + 'E':'S', + 'B':'S', + 'T':'L', + 'S':'L', + '-':'L', + + 'b':'S', + 'C':'L' + } + self.SSE_name = { + 'H' : 'Alpha helix', + 'G' : '3-10 helix', + 'I' : 'Pi helix', + + 'E' : 'Extended strand', + 'B' : 'Isolated beta-bridge', + + 'T' : 'Turn', + 'S' : 'Bend', + '-' : 'None', + + 'b' : 'Isolated beta-bridge', + 'C' : 'Coil' + } + + self.SSE_col_RGB = { + 'H':(255, 0, 0), + 'G':(255, 0,128), + 'I':(255,170,170), + + 'E':(255,255, 0), + 'B':(153,119, 85), + + 'T':( 0,255,255), + 'S':(153,255,119), + '-':(179,179,179), + + 'b':(153,119, 85), # same as 'B' + 'C':( 0, 0,255) + } + self.SSE_col = {} + for sse in self.SSE_col_RGB.keys(): + self.SSE_col[sse] = '#%s%s%s' % (hex(self.SSE_col_RGB[sse][0])[2:].zfill(2), + hex(self.SSE_col_RGB[sse][1])[2:].zfill(2), + hex(self.SSE_col_RGB[sse][2])[2:].zfill(2)) + + self.SSE_res_dict = {} + self.SSE_sel_dict = {} + + w = Tkinter.Label(self.dialog.interior(), +# text = '\nDSSP Plugin for PyMOL\nHongbo Zhu, 2011.\n\nColor proteins according to the secondary structure determined by DSSP.', + text = '\nDSSP and Stride Plugin for PyMOL\nby Hongbo Zhu, 2011\n', + background = 'black', foreground = 'green' + ) + w.pack(expand = 1, fill = 'both', padx = 10, pady = 5) + + # make a few tabs within the dialog + self.notebook = Pmw.NoteBook(self.dialog.interior()) + self.notebook.pack(fill = 'both', expand=1, padx=10, pady=10) + + + ###################### + # Tab : Structure Tab + ###################### + page = self.notebook.add('Structure') + self.notebook.tab('Structure').focus_set() + group_struc = Tkinter.LabelFrame(page, text = 'Structure') + group_struc.pack(fill='both', expand=True, padx=10, pady=5) + + pymol_sel_ent = Pmw.EntryField(group_struc, + label_text='PyMOL selection/object:', + labelpos='wn', + entry_textvariable=self.pymol_sel + ) + + dssp_bin_ent = Pmw.EntryField(group_struc, + label_text='DSSP binary:', labelpos='wn', + entry_textvariable=self.dssp_bin, + entry_width=20) + dssp_bin_but = Tkinter.Button(group_struc, text = 'Browse...', + command = self.getDSSPBin) + + stride_bin_ent = Pmw.EntryField(group_struc, + label_text='Stride binary:', labelpos='wn', + entry_textvariable=self.stride_bin, + entry_width=20) + stride_bin_but = Tkinter.Button(group_struc, text = 'Browse...', + command = self.getStrideBin) + + # arrange widgets using grid + pymol_sel_ent.grid(sticky='we', row=0, column=0, + columnspan=2, padx=5, pady=5) + dssp_bin_ent.grid(sticky='we', row=1, column=0, padx=5, pady=1) + dssp_bin_but.grid(sticky='we', row=1, column=1, padx=5, pady=1) + stride_bin_ent.grid(sticky='we', row=2, column=0, padx=5, pady=1) + stride_bin_but.grid(sticky='we', row=2, column=1, padx=5, pady=1) + group_struc.columnconfigure(0, weight=9) + group_struc.columnconfigure(1, weight=1) + + + ###################### + # Tab : Color Tab + ###################### + # H = alpha helix + # G = 3-helix (3/10 helix) + # I = 5 helix (pi helix) + + # E = extended strand, participates in beta ladder + # B = residue in isolated beta-bridge + + # T = hydrogen bonded turn + # S = bend + # + # H Alpha helix (4-12) + # G 3-10 helix + # I pi helix + # E Strand + # B Isolated beta-bridge residue + # T Turn + # S Bend + # - None + + page = self.notebook.add('Color') + + group_sse_color = Tkinter.LabelFrame(page, text = 'Secondary Structure Element Color') + group_sse_color.grid(sticky='eswn', row=0, column=0, padx=10, pady=5) + + # colors for DSSP surface + H_col_lab = Tkinter.Label(group_sse_color, text='Alpha helix (H):') + self.H_col_but = Tkinter.Button(group_sse_color, + bg=self.SSE_col['H'], + activebackground=self.SSE_col['H'], + command = self.custermizeHColor) + G_col_lab = Tkinter.Label(group_sse_color, text='3-10 helix (G):') + self.G_col_but = Tkinter.Button(group_sse_color, + bg=self.SSE_col['G'], + activebackground=self.SSE_col['G'], + command = self.custermizeGColor) + I_col_lab = Tkinter.Label(group_sse_color, text='PI helix (I):') + self.I_col_but = Tkinter.Button(group_sse_color, + bg=self.SSE_col['I'], + activebackground=self.SSE_col['G'], + command = self.custermizeIColor) + + + E_col_lab = Tkinter.Label(group_sse_color, text='Extended strand (E):') + self.E_col_but = Tkinter.Button(group_sse_color, + bg=self.SSE_col['E'], + activebackground=self.SSE_col['E'], + command = self.custermizeEColor) + B_col_lab = Tkinter.Label(group_sse_color, text='Isolated beta-bridge (B):') + self.B_col_but = Tkinter.Button(group_sse_color, + bg=self.SSE_col['B'], + activebackground=self.SSE_col['B'], + command = self.custermizeBColor) + + T_col_lab = Tkinter.Label(group_sse_color, text='Turn (T):') + self.T_col_but = Tkinter.Button(group_sse_color, + bg=self.SSE_col['T'], + activebackground=self.SSE_col['T'], + command = self.custermizeTColor) + S_col_lab = Tkinter.Label(group_sse_color, text='Bend (DSSP S):') + self.S_col_but = Tkinter.Button(group_sse_color, + bg=self.SSE_col['S'], + activebackground=self.SSE_col['S'], + command = self.custermizeSColor) + N_col_lab = Tkinter.Label(group_sse_color, text='None (DSSP NA):') + self.N_col_but = Tkinter.Button(group_sse_color, + bg=self.SSE_col['-'], + activebackground=self.SSE_col['-'], + command = self.custermizeNColor) + + + b_col_lab = Tkinter.Label(group_sse_color, text='Isolated beta-bridge (Stride b):') + self.b_col_but = Tkinter.Button(group_sse_color, + bg=self.SSE_col['b'], + activebackground=self.SSE_col['b'], + command = self.custermizebColor) + + C_col_lab = Tkinter.Label(group_sse_color, text='Coil (Stride C):') + self.C_col_but = Tkinter.Button(group_sse_color, + bg=self.SSE_col['C'], + activebackground=self.SSE_col['C'], + command = self.custermizeCColor) + + + + H_col_lab.grid(sticky='e', row=0, column=0, padx=5, pady=3) + self.H_col_but.grid(sticky='we', row=0, column=1, padx=5, pady=3) + G_col_lab.grid(sticky='e', row=1, column=0, padx=5, pady=3) + self.G_col_but.grid(sticky='we', row=1, column=1, padx=5, pady=3) + I_col_lab.grid(sticky='e', row=2, column=0, padx=5, pady=3) + self.I_col_but.grid(sticky='we', row=2, column=1, padx=5, pady=3) + + E_col_lab.grid(sticky='e', row=0, column=2, padx=5, pady=3) + self.E_col_but.grid(sticky='we', row=0, column=3, padx=5, pady=3) + B_col_lab.grid(sticky='e', row=1, column=2, padx=5, pady=3) + self.B_col_but.grid(sticky='we', row=1, column=3, padx=5, pady=3) + + T_col_lab.grid(sticky='e', row=0, column=4, padx=5, pady=3) + self.T_col_but.grid(sticky='we', row=0, column=5, padx=5, pady=3) + S_col_lab.grid(sticky='e', row=1, column=4, padx=5, pady=3) + self.S_col_but.grid(sticky='we', row=1, column=5, padx=5, pady=3) + N_col_lab.grid(sticky='e', row=2, column=4, padx=5, pady=3) + self.N_col_but.grid(sticky='we', row=2, column=5, padx=5, pady=3) + + b_col_lab.grid(sticky='e', row=2, column=2, padx=5, pady=3) + self.b_col_but.grid(sticky='we', row=2, column=3, padx=5, pady=3) + + C_col_lab.grid(sticky='e', row=3, column=4, padx=5, pady=3) + self.C_col_but.grid(sticky='we', row=3, column=5, padx=5, pady=3) + ###################### + # Tab : About Tab + ###################### + page = self.notebook.add('About') + group_about = Tkinter.LabelFrame(page, text = 'About DSSP and Stride Plugin for PyMOL') + group_about.grid(sticky='we', row=0,column=0,padx=5,pady=3) + about_plugin = """ Assign and color secondary structures using DSSP or Stride. +by Hongbo Zhu +Please cite this plugin if you use it in a publication. +Hongbo Zhu. DSSP and Stride plugin for PyMOL, 2011, BIOTEC, TU Dresden. +""" + + label_about = Tkinter.Label(group_about,text=about_plugin) + label_about.grid(sticky='we', row=0, column=0, padx=5, pady=10) + + self.notebook.setnaturalsize() + + return + + + def getDSSPBin(self): + dssp_bin_fname = tkFileDialog.askopenfilename( + title='DSSP Binary', initialdir='', + filetypes=[('all','*')], parent=self.parent) + if dssp_bin_fname: # if nonempty + self.dssp_bin.set(dssp_bin_fname) + return + + def getStrideBin(self): + stride_bin_fname = tkFileDialog.askopenfilename( + title='Stride Binary', initialdir='', + filetypes=[('all','*')], parent=self.parent) + if stride_bin_fname: # if nonempty + self.stride_bin.set(stride_bin_fname) + return + + def runDSSPOneObj(self, one_obj_sel): + """ Run DSSP on only one object. + + @param one_obj_sel: the selection/object involving only one object + @param type: string + """ + SSE_res = { + 'H':{}, 'G':{}, 'I':{}, + 'E':{}, 'B':{}, + 'T':{}, 'S':{}, '-':{} + } + SSE_sel = { + 'H':None, 'G':None, 'I':None, + 'E':None, 'B':None, + 'T':None, 'S':None, '-':None + } + + pdb_fn = None + pdb_os_fh, pdb_fn = tempfile.mkstemp(suffix='.pdb') # file os handle, file name + os.close(pdb_os_fh) + # DSSP 2.0.4 ignores all residues after 1st TER in the same chain + v = cmd.get(name='pdb_use_ter_records') + if v: cmd.set(name='pdb_use_ter_records', value=0) # do not insert TER into the pdb + cmd.save(filename=pdb_fn, selection=one_obj_sel) + if v: cmd.set(name='pdb_use_ter_records', value=v) # restore old value + + if VERBOSE: + print 'Selection %s saved to %s.' % (one_obj_sel, pdb_fn) + + if pdb_fn is None: + print 'WARNING: DSSP has no pdb file to work on!' + return None + + print 'Running DSSP for %s ...' % (one_obj_sel,) + dssp_sse_dict = {} + if sys.version_info >= (2,4): + dssp_proc = subprocess.Popen([self.dssp_bin.get(), pdb_fn], + stdout=subprocess.PIPE, + stderr=subprocess.STDOUT) + dssp_stdout, dssp_stderr = dssp_proc.communicate() + else: # use os.system + tempfile + dssp_tmpout_os_fh, dssp_tmpout_fn = tempfile.mkstemp(suffix='.dssp') + os.close(dssp_tmpout_os_fh) + dssp_cmd = '%s %s > %s' % (self.dssp_bin.get(), pdb_fn, dssp_tmpout_fn) + os.system(dssp_cmd) + fh = open(dssp_tmpout_fn) + dssp_stdout = ''.join(fh.readlines()) + fh.close() + + sse_started = False + for line in dssp_stdout.splitlines(): + if line.startswith(' # RESIDUE'): + sse_started = True + continue + elif line.startswith(' !!!'): + sse_started = False + continue + elif sse_started: + if len(line) < 10 or line[9] == ' ': continue + ch,resname = line[11],line[13] + residen,sscode = line[5:11].strip(),line[16] # residen = resnum+icode, col 10 is for icode + if sscode == ' ': sscode = '-' + k = (ch,resname,residen) + dssp_sse_dict[k] = sscode + + self.dssp_rlt_dict[one_obj_sel] = dssp_sse_dict + print 'Got SSE for %d residues.' % (len(self.dssp_rlt_dict[one_obj_sel]),) + + # group residues according to their SSE, and chain name + for k in self.dssp_rlt_dict[one_obj_sel].keys(): + #res = '/%s//%s/%d%s/' % (sel,k[0],k[1][1], k[1][2].strip()) # sel name, chain ID, res serial num, icode + sse = self.dssp_rlt_dict[one_obj_sel][k] + chn = k[0] + res = k[2] + if VERBOSE: print '(%s) and \"%s\"/%s/ sse=%s' % (str(one_obj_sel), + chn.strip(), res, sse) + SSE_res[sse].setdefault(chn,[]).append(res) + + self.SSE_res_dict[one_obj_sel] = SSE_res + + for sse in self.DSSP_SSE_list: + sse_sel_name = self.selectSSE(one_obj_sel, sse) + SSE_sel[sse] = sse_sel_name + + self.SSE_sel_dict[one_obj_sel] = SSE_sel + + print '\nNumber of residues with SSE element:' + for sse in self.DSSP_SSE_list: + num = sum([len(SSE_res[sse][chn]) for chn in SSE_res[sse]]) + print '%20s (%s) : %5d' % (self.SSE_name[sse], sse, num) + print + + # clean up pdb_fn and dssp_tmpout_fn created by tempfile.mkstemp() + if os.path.isfile(pdb_fn): + os.remove(pdb_fn) + if sys.version_info < (2,4) and os.path.isfile(dssp_tmpout_fn): + os.remove(dssp_tmpout_fn) + + return + + + def runDSSP(self): + """ + @return: whether DSSP has been executed successfully + @rtype: boolean + """ + # delete old results + self.sel_obj_list = [] + self.dssp_rlt_dict = {} + self.SSE_res_dict = {} + self.SSE_sel_dict = {} + + pdb_fn = None + sel_name= None + sel = self.pymol_sel.get() + + if len(sel) > 0: # if any pymol selection/object is specified + # save the pymol selection/object in the tmp dir + all_sel_names = cmd.get_names('all') # get names of all selections + if sel in all_sel_names: + if cmd.count_atoms(sel) == 0: + err_msg = 'ERROR: The selection %s is empty.' % (sel,) + print 'ERROR: %s' % (err_msg,) + tkMessageBox.showinfo(title='ERROR', message=err_msg) + return False + else: + sel_name = sel + # no selection/object with the input name is found + # we assume either a single-word selector or + # some other selection-expression is used + # NOTE: if more than one selection is selected, they should be + # saved separately and SSE should be calculated for each + # of the selections. + else: + print 'The selection/object you specified is not found.' + print 'Your input will be interpreted as a selection-expression.' + # check whether the selection is empty + tmpsel = self.randomSeleName(prefix='your_sele_') + cmd.select(tmpsel,sel) + if cmd.count_atoms(tmpsel) == 0: + cmd.delete(tmpsel) + err_msg = 'ERROR: The selection %s is empty.' % (sel,) + print 'ERROR: %s' % (err_msg,) + tkMessageBox.showinfo(title='ERROR', message=err_msg) + return False + else: + sel_name = tmpsel + else: # what structure do you want DSSP to work on? + err_msg = 'No PyMOL selection/object specified!' + print 'ERROR: %s' % (err_msg,) + tkMessageBox.showinfo(title='ERROR', message=err_msg) + return False + + # each object in the selection is treated as an independent struc + objlist = cmd.get_object_list(sel_name) + self.ss_asgn_prog = 'DSSP' + print 'Starting %s ...' % (self.ss_asgn_prog, ) + + for objname in objlist: + self.sel_obj_list.append('%s and %s' % (sel_name, objname)) + self.runDSSPOneObj(self.sel_obj_list[-1]) + +## cmd.delete(tmpsel) + + return True + + + def runStrideOneObj(self, one_obj_sel): + """ Run Stride on only one object. + + @param one_obj_sel: the selection/object involving only one object + @param type: string + """ + SSE_res = { + 'H':{}, 'G':{}, 'I':{}, + 'E':{}, 'B':{}, 'b':{}, + 'T':{}, 'C':{} + } + SSE_sel = { + 'H':None, 'G':None, 'I':None, + 'E':None, 'B':None, 'b':None, + 'T':None, 'C':None + } + + pdb_fn = None + pdb_os_fh, pdb_fn = tempfile.mkstemp(suffix='.pdb') # file os handle, file name + os.close(pdb_os_fh) + cmd.save(filename=pdb_fn, selection=one_obj_sel) + if VERBOSE: + print 'Selection %s saved to %s.' % (one_obj_sel, pdb_fn) + + if pdb_fn is None: + print 'WARNING: Stride has no pdb file to work on!' + return None + + print 'Running Stride for %s ...' % (one_obj_sel,) + stride_sse_dict = {} + if sys.version_info >= (2,4): + stride_proc = subprocess.Popen([self.stride_bin.get(), pdb_fn], + stdout=subprocess.PIPE, + stderr=subprocess.STDOUT + ) + stride_stdout, stride_stderr = stride_proc.communicate() + else: # use os.system + tempfile + stride_tmpout_os_fh, stride_tmpout_fn = tempfile.mkstemp(suffix='.stride') + os.close(stride_tmpout_os_fh) + stride_cmd = '%s %s > %s' % (self.stride_bin.get(), pdb_fn, stride_tmpout_fn) + os.system(stride_cmd) + fh = open(stride_tmpout_fn) + stride_stdout = ''.join(fh.readlines()) + fh.close() + + for line in stride_stdout.splitlines(): + if line.startswith('ASG'): + resname,ch=line[5:8], line[9] + # according to stride doc, col 12-15 are used for residue number + # actually, resnum+icode is used in 12-15. + # In addition, if residue number is 4-digit or longer, + # the field 12-16 or more are used. + + if line[15] == ' ': # col 16 is not occupied + residen, sscode = line[11:15].strip(),line[24] # residen = resnum+icode + else: + shift = line[15:].find(' ') # find where residen stops + residen, sscode = line[11:15+shift].strip(),line[24+shift] + + k = (ch,resname,residen) + stride_sse_dict[k] = sscode + + self.stride_rlt_dict[one_obj_sel] = stride_sse_dict + print 'Got SSE for %d residues.' % (len(self.stride_rlt_dict[one_obj_sel]),) + + # group residues according to their SSE, and chain name + for k in self.stride_rlt_dict[one_obj_sel].keys(): + sse=self.stride_rlt_dict[one_obj_sel][k] + chn = k[0] # this can be a space! + res = k[2] + if VERBOSE: print '(%s) and \"%s\"/%s/ sse=%s' % (str(one_obj_sel), + chn.strip(), res, sse) + SSE_res[sse].setdefault(chn,[]).append(res) + + self.SSE_res_dict[one_obj_sel] = SSE_res + + for sse in self.STRIDE_SSE_list: + sse_sel_name = self.selectSSE(one_obj_sel, sse) + SSE_sel[sse] = sse_sel_name + self.SSE_sel_dict[one_obj_sel] = SSE_sel + + print '\nNumber of residues with SSE element:' + for sse in self.STRIDE_SSE_list: + num = sum([len(SSE_res[sse][chn]) for chn in SSE_res[sse]]) + print '%20s (%s) : %5d' % (self.SSE_name[sse], sse, num) + print + + # clean up pdb_fn and dssp_tmpout_fn created by tempfile.mkstemp() + if os.path.isfile(pdb_fn): + os.remove(pdb_fn) + if sys.version_info < (2,4) and os.path.isfile(stride_tmpout_fn): + os.remove(stride_tmpout_fn) + + return True + + + def runStride(self): + """ + """ + # delete old results + self.sel_obj_list = [] + self.stride_rlt_dict = {} + self.SSE_res_dict = {} + self.SSE_sel_dict = {} + + pdb_fn = None + sel_name= None + sel = self.pymol_sel.get() + + if len(sel) > 0: # if any pymol selection/object is specified + all_sel_names = cmd.get_names('all') # get names of all selections + if sel in all_sel_names: + if cmd.count_atoms(sel) == 0: + err_msg = 'ERROR: The selection %s is empty.' % (sel,) + print 'ERROR: %s' % (err_msg,) + tkMessageBox.showinfo(title='ERROR', message=err_msg) + return False + else: + sel_name = sel + # no selection/object with the input name is found + # we assume either a single-word selector or + # some other selection-expression is uesd + else: + print 'The selection/object you specified is not found.' + print 'Your input will be interpreted as a selection-expression.' + tmpsel = self.randomSeleName(prefix='your_sele_') + cmd.select(tmpsel,sel) + if cmd.count_atoms(tmpsel) == 0: + cmd.delete(tmpsel) + err_msg = 'ERROR: The selection %s is empty.' % (sel,) + print 'ERROR: %s' % (err_msg,) + tkMessageBox.showinfo(title='ERROR', message=err_msg) + return False + else: + sel_name = tmpsel + + else: # what structure do you want Stride to work on? + err_msg = 'No PyMOL selection/object specified!' + print 'ERROR: %s' % (err_msg,) + tkMessageBox.showinfo(title='ERROR', message=err_msg) + return False + + # each object in the selection is treated as an independent struc + objlist = cmd.get_object_list(sel_name) + self.ss_asgn_prog = 'Stride' + print 'Starting %s ...' % (self.ss_asgn_prog, ) + + for objname in objlist: + self.sel_obj_list.append('%s and %s' % (sel_name, objname)) + self.runStrideOneObj(self.sel_obj_list[-1]) + + return True + + + def randomSeleName(self, prefix='sele',suffix=''): + """ generate a random selection name. + """ + sel_list = cmd.get_names('all') + sel_dict = dict(zip(sel_list, range(len(sel_list)))) + sel_name = '%s%d%s' % (prefix, random.randint(1000,9999), suffix) + while(sel_name in sel_dict): + sel_name = '%s%d%s' % (prefix, random.randint(1000,9999), suffix) + + return sel_name + + + def selectSSE(self, sel, sse): + """ generate selector for selecting all residues having the given sse. + return the selection name. + """ + #sel = self.pymol_sel.get() + sel_list_chn = [] + + if VERBOSE: print '\nSelecting SSE %s ... \n' % (sse) + + for chn in self.SSE_res_dict[sel][sse]: # color one chain at a time + if chn == ' ': chn_str = '-' + else: chn_str = chn + if VERBOSE: print 'Selecting SSE %s on chain %s ... \n' % (sse, chn) + limit=150 # color every 150 residues + sel_name_chn = self.randomSeleName(prefix='%s_%s_%s_' % ('_'.join(sel.split()),chn_str,sse)) + + if len(self.SSE_res_dict[sel][sse][chn]) < limit: + #sel_expr = '/%s//%s/%s/' % (sel,chn.strip(), '+'.join(self.SSE_res[sse][chn])) + # always quote chain name in case it is empty (otherwise sel misinterpreted) + sel_expr = '(%s) and \"%s\"/%s/' % (sel,chn.strip(), '+'.join(self.SSE_res_dict[sel][sse][chn])) + cmd.select(sel_name_chn, sel_expr) + if VERBOSE: print 'select %s, %s' % (sel_name_chn, sel_expr) + sel_list_chn.append(sel_name_chn) + else: + rn = len(self.SSE_res_dict[sel][sse][chn]) + print 'total number of res with %s = %d' % (sse,rn) + sz = int(math.ceil(rn/float(limit))) + sel_list_seg = [] + for i in xrange(sz): + s,e = i*limit, min((i+1)*limit, rn) + print s,e + #sel_expr = '/%s//%s/%s/' % (sel,chn.strip(), '+'.join(self.SSE_res[sse][chn][s:e])) + # always quote chain name in case it is empty (otherwise sel misinterpreted) + sel_expr = '(%s) and \"%s\"/%s/' % (sel, chn.strip(), + '+'.join(self.SSE_res_dict[sel][sse][chn][s:e])) + sel_name_seg = self.randomSeleName(prefix='%s_%s_%s_tmp_' % ('_'.join(sel.split()),chn_str,sse)) + cmd.select(sel_name_seg, sel_expr) + if VERBOSE: print 'select %s, %s' % (sel_name_seg, sel_expr) + sel_list_seg.append(sel_name_seg) + + sel_expr = ' or '.join(sel_list_seg) + cmd.select(sel_name_chn, sel_expr) + if VERBOSE: print 'select %s, %s' % (sel_name_chn, sel_expr) + [cmd.delete(asel) for asel in sel_list_seg] + sel_list_chn.append(sel_name_chn) + + if len(sel_list_chn) > 0: + sel_name = self.randomSeleName(prefix='%s_%s_%s_' % ('_'.join(sel.split()), sse, self.ss_asgn_prog)) + sel_expr = ' or '.join(sel_list_chn) + cmd.select(sel_name,sel_expr) + [cmd.delete(asel) for asel in sel_list_chn] + else: + print 'INFO: No residues are assigned to SSE \'%s\'.' % (sse,) + sel_name = None + + return sel_name + + + def updateColor(self): + if self.ss_asgn_prog is None: + err_msg = 'Run DSSP or Stride to assign secondary structures first!' + print 'ERROR: %s' % (err_msg,) + tkMessageBox.showinfo(title='ERROR', message=err_msg) + else: + print 'Update color for %s' % (self.pymol_sel.get()), + print 'using secondary structure assignment by %s' % (self.ss_asgn_prog,) + + if self.ss_asgn_prog == 'DSSP': SSE_list = self.DSSP_SSE_list + elif self.ss_asgn_prog == 'Stride': SSE_list = self.STRIDE_SSE_list + + for sse in SSE_list: # give color names + cmd.set_color('%s_color'%(sse,), self.SSE_col_RGB[sse]) + for sse in SSE_list: # color each SSE + for sel_obj in self.sel_obj_list: + if self.SSE_sel_dict[sel_obj][sse] is not None: + cmd.color('%s_color'%(sse,), self.SSE_sel_dict[sel_obj][sse]) + print 'color',self.SSE_sel_dict[sel_obj][sse],',',self.SSE_col_RGB[sse] + else: + print 'No residues with SSE \'%s\' to color.' % (sse,) + + return + + def updateSS(self): + if self.ss_asgn_prog is None: + err_msg = 'Run DSSP or Stride to assign secondary structures first!' + print 'ERROR: %s' % (err_msg,) + tkMessageBox.showinfo(title='ERROR', message=err_msg) + else: + print 'Update secondary structures for %s' % (self.pymol_sel.get()), + print 'using secondary structure assignment by %s' % (self.ss_asgn_prog,) + print 'SSE mapping: (H,G,I) ==> H; (E,B,b) ==> S; (T,S,-,C) ==> L' + + if self.ss_asgn_prog == 'DSSP': SSE_list = self.DSSP_SSE_list + elif self.ss_asgn_prog == 'Stride': SSE_list = self.STRIDE_SSE_list + + for sse in SSE_list: + for sel_obj in self.sel_obj_list: + if self.SSE_sel_dict[sel_obj][sse] is not None: + cmd.alter(self.SSE_sel_dict[sel_obj][sse], 'ss=\'%s\''% (self.SSE_map[sse],)) + print 'alter %s, ss=%s' % (self.SSE_sel_dict[sel_obj][sse], self.SSE_map[sse]) + else: + print 'No residue with SSE %s to update ss.' % (sse,) + + # show cartoon for the input selection, and rebuild + cmd.show('cartoon',self.pymol_sel.get()) + cmd.rebuild(self.pymol_sel.get()) + return + + + def custermizeHColor(self): + self.custermizeSSEColor('H') + def custermizeGColor(self): + self.custermizeSSEColor('G') + def custermizeIColor(self): + self.custermizeSSEColor('I') + + def custermizeEColor(self): + self.custermizeSSEColor('E') + def custermizeBColor(self): + self.custermizeSSEColor('B') + + def custermizeTColor(self): + self.custermizeSSEColor('T') + def custermizeSColor(self): + self.custermizeSSEColor('S') + def custermizeNColor(self): + self.custermizeSSEColor('-') + + def custermizebColor(self): + self.custermizeSSEColor('b') + def custermizeCColor(self): + self.custermizeSSEColor('C') + + def custermizeSSEColor(self,sse): + SSE_col_but = { + 'H':self.H_col_but, + 'G':self.G_col_but, + 'I':self.I_col_but, + + 'E':self.E_col_but, + 'B':self.B_col_but, + + 'T':self.T_col_but, + 'S':self.S_col_but, + '-':self.N_col_but, + + 'b':self.b_col_but, + 'C':self.C_col_but, + } + try: + color_tuple, color = tkColorChooser.askcolor(color=self.SSE_col[sse]) + if color_tuple is not None and color is not None: + self.SSE_col_RGB[sse] = color_tuple + self.SSE_col[sse] = color + + SSE_col_but[sse]['bg']=self.SSE_col[sse] + SSE_col_but[sse]['activebackground']=self.SSE_col[sse] + SSE_col_but[sse].update() + except Tkinter._tkinter.TclError: + print 'Old color (%s) will be used.' % (self.mesh_col) + + + def execute(self, butcmd): + """ Run the cmd represented by the botton clicked by user. + """ + if butcmd == 'OK': + print 'is everything OK?' + + elif butcmd == 'Run DSSP': + rtn = self.runDSSP() + if rtn and VERBOSE: print 'Done with DSSP!' + + elif butcmd == 'Run Stride': + rtn = self.runStride() + if rtn and VERBOSE: print 'Done with Stride!' + + elif butcmd == 'Update Color': + self.updateColor() + + elif butcmd == 'Update ss': + self.updateSS() + + elif butcmd == 'Exit': + print 'Exiting DSSP and Stride Plugin ...' + if __name__ == '__main__': + self.parent.destroy() + else: + self.dialog.withdraw() + print 'Done.' + else: + print 'Exiting DSSP and Stride Plugin because of unknown button click ...' + self.dialog.withdraw() + print 'Done.' + + + def quit(self): + self.dialog.destroy() + + +############################################# +# +# +# Create demo in root window for testing. +# +# +############################################## +if __name__ == '__main__': + + class App: + def my_show(self,*args,**kwargs): + pass + + app = App() + app.root = Tkinter.Tk() + Pmw.initialise(app.root) + app.root.title('It seems to work!') + + widget = DSSPPlugin(app) + app.root.mainloop() diff --git a/.pymol/startup/dssp_stride.pyc b/.pymol/startup/dssp_stride.pyc new file mode 100644 index 0000000..75c888b Binary files /dev/null and b/.pymol/startup/dssp_stride.pyc differ diff --git a/.pymol/startup/propka.py b/.pymol/startup/propka.py new file mode 100644 index 0000000..a8c4d0b --- /dev/null +++ b/.pymol/startup/propka.py @@ -0,0 +1,887 @@ +''' +Described at PyMOL wiki: +http://www.pymolwiki.org/index.php/propka + +#------------------------------------------------------------------------------- +# Name: propka for pymol +# Purpose: To fetch and display the pka values for protein of intetest +# +# Author: Troels E. Linnet +# +# Created: 14/08/2011 +# Copyright: (c) Troels E. Linnet 2011 +# Contact: tlinnet snabela gmail dot com +# Licence: Free for all +# +# +#------------------------------------------------------------------------------- + + The PROPKA method is developed by the + Jensen Research Group + Department of Chemistry + University of Copenhagen + + Please cite these references in publications: +Hui Li, Andrew D. Robertson, and Jan H. Jensen +"Very Fast Empirical Prediction and Interpretation of Protein pKa Values" +Proteins, 2005, 61, 704-721. + +Delphine C. Bas, David M. Rogers, and Jan H. Jensen +"Very Fast Prediction and Rationalization of pKa Values for Protein-Ligand Complexes" +Proteins, 2008, 73, 765-783. + +Mats H.M. Olsson, Chresten R. Soendergard, Michal Rostkowski, and Jan H. Jensen +"PROPKA3: Consistent Treatment of Internal and Surface Residues in Empirical pKa predictions" +Journal of Chemical Theory and Computation, 2011 7 (2), 525-537 + +Chresten R. Soendergaard, Mats H.M. Olsson, Michaz 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, 2011 in press +""" +#------------------------------------------------------------------------------- +# The script needs mechanize to run. +# On windows, it is not easy to make additional modules available for pymol. So put in into your working folder. +#1)The easy manual way: +#a)Go to: http://wwwsearch.sourceforge.net/mechanize/download.html +#b)Download mechanize-0.2.5.zip. http://pypi.python.org/packages/source/m/mechanize/mechanize-0.2.5.zip +#c)Extract to .\mechanize-0.2.5 then move the in-side folder "mechanize" to your folder with propka.py. The rest of .\mechanize-0.2.5 you don't need. +#You can also see other places where you could put the "mechanize" folder. Write this in pymol to see the paths where pymol is searching for "mechanize" +# import sys; print(sys.path) + +#------------------------------------------------------------------------------- +""" +Example for pymol script to start the functions. For example: trypropka.pml +Execute with pymol or start pymol and: File->Run->trypropka.pml +############################################################################################################################################################################################################################## + +### Point to your directory with your pdb file and where to save the results +#cd /homes/linnet/Documents/Speciale/5NT-project/Mutant-construct/predict_reactivity/propka +cd C:/Users/tlinnet/Documents/My Dropbox/Speciale/5NT-project/Mutant-construct/predict_reactivity/propka + +### The fastest method is just to write propka. Then the last pymol molecule is assumed and send to server. verbose=yes makes the script gossip mode. +import propka + +fetch 4ins, async=0 +propka +#fetch 1hp1, async=0 +#propka logtime=_, resi=5-10.20-30, resn=CYS.ATP.TRP, verbose=yes + +### Fetch 4ins from web. async make sure, we dont execute script before molecule is loaded. The resi and resn prints the interesting results right to command line. +#fetch 4ins, async=0 +#propka chain=*, resi=5-10.20-30, resn=ASP.CYS, logtime=_ + +### If there is no web connection, one can process a local .pka file. Either from a previous run or from a downloaded propka webpage result. +### Then run and point to .pka file with: pkafile=./Results_propka/pkafile.pka Remember the dot "." in the start, to make it start in the current directory. +#load 4ins.pdb +#propka pkafile=./Results_propka/4ins_.pka, resi=18.25-30, resn=cys, + +### Some more examples. This molecule has 550 residues, so takes a longer time. We select to run the last molecule, by writing: molecule=1hp1 +#fetch 4ins, async=0 +#fetch 1hp1, async=0 +#propka molecule=1hp1, chain=A, resi=300-308.513, resn=CYS.ATP.TRP, logtime=_, verbose=no, showresult=no +#propka molecule=1hp1, pkafile=./Results_propka/1hp1_.pka, verbose=yes + +### One can also just make a lookup for a protein. Use function: getpropka +### Note. This does only print the result to the pymol command line +#getpropka source=ID, PDBID=4ake, logtime=_, showresult=yes +#getpropka source=ID, PDBID=4ins, logtime=_, server_wait=10.0, verbose=yes, showresult=no +############################################Input parameters: propka############################################ +############# The order of input and changable things: +############# propka(molecule="NIL",chain="*",resi="0",resn="NIL",method="upload",logtime=time.strftime("%m%d",time.localtime()),server_wait=3.0,version="v3.1",verbose="no",showresult="no",pkafile="NIL") +# method : method=upload is default. This sends .pdb file and request result from propka server. +## method=file will only process a manual .pka file, and write a pymol command file. No use of mechanize. +## If one points to an local .pka file, then method is auto-changed to method=file. This is handsome in off-line environment, ex. teaching or seminar. +# pkafile: Write the path to .pka file. Ex: pkafile=./Results_propka/4ins_.pka +# molecule : name of the molecule. Ending of file is assumed to be .pdb +# chain : which chains are saved to file, before molecule file is send to server. Separate with "." Ex: chain=A.b +# resi : Select by residue number, which residues should be printed to screen and saved to the log file: /Results_propka/_Results.log. +## Separate with "." or make ranges with "-". Ex: resi=35.40-50 +# resn : Select by residue name, which residues should be printed to screen and saved to the log file: /Results_propka/_Results.log. +## Separate with "." Ex: resn=cys.tyr +# logtime : Each execution give a set of files with the job id=logtime. If logtime is not provided, the current time is used. +## Normal it usefull to set it empty. Ex: logtime=_ +# verbose : Verbose is switch, to turn on messages for the mechanize section. This is handsome to see how mechanize works, and for error searching. +# showresult : Switch, to turn on all results in pymol command window. Ex: showresult=yes +# server_wait=10.0 is default. This defines how long time between asking the server for a result. Set no lower than 3 seconds. +# version=v3.1 is default. This is what version of propka which would be used. +## Possible: 'v3.1','v3.0','v2.0'. If a newer version is available than the current v3.1, a error message is raised to make user update the script. +############################################Input parameters: getpropka############################################ +############# The order of input and changable things: +############# getpropka(PDB="NIL",chain="*",resi="0",resn="NIL",source="upload",PDBID="",logtime=time.strftime("%Y%m%d%H%M%S",time.localtime()),server_wait=3.0,version="v3.1",verbose="no",showresult="no") +# PDB: points the path to a .pdb file. This is auto-set from propka function. +# source : source=upload is default and is set at the propka webpage. +# source=ID, PDBID=4ake , one can print to the command line, the pka value for any official pdb ID. No files are displayed in pymol. +# PDBID: is used as the 4 number/letter pdb code, when invoking source=ID. + +############################################################################################################################################################################################################################## +''' + +from __future__ import print_function + +try: + from pymol import cmd + runningpymol = 'yes' +except: + runningpymol = 'no' + pass +import time +import platform +import os +import sys + + +def propka(molecule="NIL", chain="*", resi="0", resn="NIL", method="upload", logtime=time.strftime("%m%d", time.localtime()), server_wait=3.0, version="v3.1", verbose="no", showresult="no", pkafile="NIL", makebonds="yes"): + Script_Version = "20110823" + # First we have to be sure, we give reasonable arguments + if pkafile != "NIL": + method = 'file' + assert method in ['upload', 'file'], "'method' has to be either: method=upload or method=file" + # If molecule="all", then try to get the last molecule + ##assert molecule not in ['NIL'], "You always have to provide molecule name. Example: molecule=4ins" + if molecule == "NIL": + assert len(cmd.get_names()) != 0, "Did you forget to load a molecule? There are no objects in pymol." + molecule = cmd.get_names()[-1] + # To print out to screen for selected residues. Can be separated with "." or make ranges with "-". Example: resi="4-8.10" + if resi != "0": + resi_range = ResiRange(resi) + else: + resi_range = [] + # Also works for residue names. They are all converted to bigger letters. Example: resn="cys.Tyr" + if resn != "NIL": + resn_range = ResnRange(resn) + else: + resn_range = resn + # Make chain range, and upper case. + chain = ChainRange(chain) + # Make result directory. We also the absolut path to the new directory. + Newdir = createdirs() + if method == "upload": + # We try to load mechanize. If this fail, one can always get the .pka file manual and the run: method=file + try: + from modules import mechanize + importedmechanize = 'yes' + except ImportError: + print("Import error. Is a module missing?") + print(sys.exc_info()) + print("Look if missing module is in your python path\n%s" % sys.path) + importedmechanize = 'no' + import modules.mechanize as mechanize + # The name for the new molecule + newmolecule = "%s%s" % (molecule, logtime) + # Create the new molecule from original loaded and for the specified chains. Save it, and disable the old molecule. + cmd.create("%s" % newmolecule, "%s and chain %s" % (molecule, chain)) + cmd.save("%s%s.pdb" % (Newdir, newmolecule), "%s" % newmolecule) + cmd.disable("%s" % molecule) + if molecule == "all": + cmd.enable("%s" % molecule) + cmd.show("cartoon", "%s" % molecule) + # Let the new molecule be shown in cartoon. + cmd.hide("everything", "%s" % newmolecule) + cmd.show("cartoon", "%s" % newmolecule) + # Make the absolut path to the newly created .pdb file. + PDB = "%s%s.pdb" % (Newdir, newmolecule) + source = "upload" + PDBID = "" + # Request server, and get the absolut path to the result file. + pkafile = getpropka(PDB, chain, resi, resn, source, PDBID, logtime, server_wait, version, verbose, showresult) + # Open the result file and put in into a handy list. + list_results, ligands_results = importpropkaresult(pkafile) + if method == "file": + assert pkafile not in ['NIL'], "You have to provide path to file. Example: pkafile=./Results_propka/4ins_2011.pka" + assert ".pka" in pkafile, 'The propka result file should end with ".pka" \nExample: pkafile=./Results_propka/4ins_2011.pka \npkafile=%s' % (pkafile) + # The name for the molecule we pass to the writing script of pymol commands + newmolecule = "%s" % molecule + cmd.hide("everything", "%s" % newmolecule) + cmd.show("cartoon", "%s" % newmolecule) + # We open the result file we have got in the manual way and put in into a handy list. + list_results, ligands_results = importpropkaresult(pkafile) + # Then we print the interesting residues to the screen. + printpropkaresult(list_results, resi, resi_range, resn, resn_range, showresult, ligands_results) + # Now create the pymol command file. This should label the protein. We get back the absolut path to the file, so we can execute it. + result_pka_pymol_name = writepymolcmd(newmolecule, pkafile, verbose, makebonds) + # Now run our command file. But only if we are running pymol. + if runningpymol == 'yes': + cmd.do("run %s" % result_pka_pymol_name) + ##if runningpymol=='yes': cmd.do("@%s"%result_pka_pymol_name) + return(list_results) +if runningpymol != 'no': + cmd.extend("propka", propka) + + +def getpropka(PDB="NIL", chain="*", resi="0", resn="NIL", source="upload", PDBID="", logtime=time.strftime("%Y%m%d%H%M%S", time.localtime()), server_wait=3.0, version="v3.1", verbose="no", showresult="no"): + try: + import modules.mechanize as mechanize + importedmechanize = 'yes' + except ImportError: + print("Import error. Is a module missing?") + print(sys.exc_info()) + print("Look if missing module is in your python path \n %s" % sys.path) + importedmechanize = 'no' + propka_v_201108 = 3.1 + url = "http://propka.ki.ku.dk/" + assert version in ['v2.0', 'v3.0', 'v3.1'], "'version' has to be either: 'v2.0', 'v3.0', 'v3.1'" + assert source in ['ID', 'upload', 'addr', 'input_file'], "'source' has to be either: 'ID', 'upload', 'addr', 'input_file'" + if source == "upload": + assert PDB not in ['NIL'], "You always have to provide PDB path. Example: PDB=.\Results_propka\4ins2011.pdb" + if source == "ID": + assert len(PDBID) == 4, "PDBID has to be 4 characters" + # To print out to screen for selected residues. Can be separated with "." or make ranges with "-". Example: resi="4-8.10" + if resi != "0": + resi_range = ResiRange(resi) + else: + resi_range = [] + # Also works for residue names. They are all converted to bigger letters. Example: resn="cys.Tyr" + if resn != "NIL": + resn_range = ResnRange(resn) + else: + resn_range = resn + # Start the browser + br = mechanize.Browser() + # We pass to the server, that we are not a browser, but this python script. Can be used for statistics at the propka server. + br.addheaders = [('User-agent', 'pythonMechanizeClient')] + # To turn on debugging messages + # br.set_debug_http(True) + # To open the start page. + page_start = br.open(url) + read_start = page_start.read() + if verbose == 'yes': + print(br.title()) + print(br.geturl()) + # To get available forms + page_forms = [f.name for f in br.forms()] + if verbose == 'yes': + print(page_forms) + # Select first form + br.select_form(name=page_forms[0]) + # Print the current selected form, so we see that we values we start with. + if verbose == 'yes': + print(br.form) + # Print the parameters of the 'version' RadioControl button and current value + if verbose == 'yes': + print(br.find_control(name='version'), br.find_control(name='version').value) + # This is to check, that the current script is "up-to-date". + propka_v_present = float(br.find_control(name='version').value[0].replace('v', '')) + if propka_v_present > propka_v_201108: + raise UserWarning('\nNew version of propka exist.\nCheck/Update your script.\nPresent:v%s > Script:v%s' % (propka_v_present, propka_v_201108)) + # Change the parameters of the 'version' radio button and then reprint the new value. Input has to be in a list [input]. + br.form['version'] = [version] + if verbose == 'yes': + print(br.find_control(name='version').value) + # Print the parameters of the 'source' RadioControl button and current value + if verbose == 'yes': + print(br.find_control(name='source'), br.find_control(name='source').value) + # Change the parameters of the 'source' radio button and then reprint the new value. Input has to be in a list [input]. + br.form['source'] = [source] + if verbose == 'yes': + print(br.find_control(name='source').value) + # This step was the must strange and took a long time. For finding the information and the double negative way. + # One have to enable the pdb button. Read more here: http://wwwsearch.sourceforge.net/old/ClientForm/ ("# All Controls may be disabled.....) + PDBID_control = br.find_control("PDBID") + PDB_control = br.find_control("PDB") + if verbose == 'yes': + print(PDBID_control.disabled, PDB_control.disabled) + if source == "ID": + PDBID_control.disabled = False + PDB_control.disabled = True + if source == "upload": + PDBID_control.disabled = True + PDB_control.disabled = False + if verbose == 'yes': + print(PDBID_control.disabled, PDB_control.disabled) + # We create the result dir, and take with us the 'path' to the result dir. + Newdir = createdirs() + # Open all the files, and assign them. + if source == "upload": + filename = PDB + if source == "ID": + filename = PDBID + files = openfiles(Newdir, filename, logtime, source) + result_pka_file = files[0] + result_input_pka_file = files[1] + result_log = files[2] + filepath = files[3] + result_pka_pkafile = files[4] + result_pka_file_stripped = files[5] + result_pka_file_bonds = files[6] + # Print the parameters of the 'PDBID' TextControl button and current value + if source == "ID" and verbose == 'yes': + print(br.find_control(name='PDBID')) + print(br.find_control(name='PDBID').value) + # Change the parameters of the 'PDBID' TextControl and then reprint the new value. Input has just to be a string. + if source == "ID": + br.form["PDBID"] = PDBID + if source == "ID" and verbose == 'yes': + print(br.find_control(name='PDBID').value) + # Print the parameters of the 'PDB' TextControl button and current value + if source == "upload" and verbose == 'yes': + print(br.find_control(name='PDB')) + print(br.find_control(name='PDB').value) + # Change the parameters of the 'PDB' FileControl and then reprint the new value. Input has just to be a string. + if source == "upload": + PDBfilename = PDB + PDBfilenamepath = PDB + if source == "upload": + br.form.add_file(open(PDBfilename), 'text/plain', PDBfilenamepath, name='PDB') + if source == "upload" and verbose == 'yes': + print(br.find_control(name='PDB')) + print(br.find_control(name='PDB').value) + # Now reprint the current selected form, so we see that we have the right values. + if verbose == 'yes': + print(br.form) + # Make "how" we would like the next request. We would like to "Click the submit button", but we have not opened the request yet. + req = br.click(type="submit", nr=0) + # Have to pass by a mechanize exception. Thats the reason for the why True + # The error was due to: br.open(req) + # mechanize._response.httperror_seek_wrapper: HTTP Error refresh: The HTTP server returned a redirect error that would lead to an infinite loop. + # The last 30x error message was: + # OK + # I haven't been able to find the refresh problem or extend the time. So we make a pass on the raised exception. + try: + print("Now sending request to server") + br.open(req) + # If there is raised an exception, we jump through to the result page after some sleep. + except mechanize.HTTPError: + # We can extract the jobid from the current browser url. + jobid = br.geturl()[32:-5] + # We notice how the script at the server presents the final result page. + url_result = url + "pka/" + jobid + ".html" + # Now we continue to try to find the result page, until we have succes. If page doesn't exist, we wait a little. + while True: + print("Result still not there. Waiting %s seconds more" % server_wait) + time.sleep(float(server_wait)) + # To pass the "break" after the exception, we make a hack, wait and then go to the result page, which is the jobid. + try: + page_result = br.open(url_result) + read_result = page_result.read() + # If we don't receive a error in getting the result page, we break out of the while loop. + break + # If the page doesn't exist yet. We go back in the while loop. + except mechanize.HTTPError: + # Wait another round + pass + # If we get a timeout, we also wait. + except mechanize.URLError: + # Wait another round + pass + htmlresult = "The detailed result is now available at: %s" % br.geturl() + print(htmlresult) + read_result = br.response().read() + # Now save the available links from the current page. But only links that satisfy the expression. + links_result = [] + for l in br.links(url_regex='http://propka.ki.ku.dk/pka'): + links_result.append(l) + # We also extract the information for neighbour bons. This is given in the url links. + bonds = [] + for l in br.links(url_regex='http://propka.ki.ku.dk/view/new_view.cgi'): + l_split = str(l).split() + lresn = l_split[2] + lresi = l_split[3] + lchain = l_split[4] + lurl = l_split[1] + lurl_split = lurl.split("&") + lresn2 = lurl_split[1] + lchain2 = lurl_split[2] + lpka = lurl_split[3] + ldesolvation = lurl_split[4] + lneighbours = lurl_split[5:] + for i in range(len(lneighbours)): + bonds.append([lresn, lresi, lchain, lresn2, lchain2, lpka, ldesolvation, lneighbours[i]]) + # Now follow the link to the .propka_input resultpage + if len(links_result) > 1: + br.follow_link(links_result[1]) + # Now get the page text for the current link + if len(links_result) > 1: + read_result1 = br.response().read() + # Save the result + if len(links_result) > 1: + result_input_pka_file.write(read_result1) + # Now follow the link to the .pka resultpage + if len(links_result) > 1: + br.back() + result_input_pka_file.close() + # Now follow first link. "Should be" available for all versions of propka. + br.follow_link(links_result[0]) + # Now get the page for the current link + read_result0 = br.response().read() + # Save the result and close file. + result_pka_file.write(read_result0) + result_pka_file.close() + # Now get the result in a list, which is sorted + list_results, ligands_results = importpropkaresult(result_pka_pkafile) + # Print to log file + result_log.write("# executed: %s \n# logtime: %s \n# source=%s \n# PDB=%s \n# chain=%s \n# PDBID=%s \n# server_wait=%s version=%s verbose=%s showresult=%s \n# resi=%s resn=%s\n# %s \n" % (time.strftime("%Y%m%d%H%M%S", time.localtime()), logtime, source, PDB, chain, PDBID, server_wait, version, verbose, showresult, resi, resn, htmlresult)) + # Print to screen + printpropkaresult(list_results, resi, resi_range, resn, resn_range, showresult, ligands_results) + # Now write to log and the stripped file + for l in list_results: + if resi != "0" and int(l[1]) in resi_range: + result_log.write("%3s %3s %s %6s %3s %5s %3s %4s %s" % (l[0], l[1], l[2], l[3], l[4], l[5], l[6], l[7], l[8]) + '\n') + if resn != "NIL" and l[0] in resn_range and int(l[1]) not in resi_range: + result_log.write("%3s %3s %s %6s %3s %5s %3s %4s %s" % (l[0], l[1], l[2], l[3], l[4], l[5], l[6], l[7], l[8]) + '\n') + result_pka_file_stripped.write("%3s %3s %s %6s %3s %5s %3s %4s %s" % (l[0], l[1], l[2], l[3], l[4], l[5], l[6], l[7], l[8]) + '\n') + for l in ligands_results: + if resn != "NIL" and l[0] in resn_range: + result_log.write("%3s %3s %s %6s %3s %5s %3s %4s %s" % (l[0], l[1], l[2], l[3], l[4], l[5], l[6], l[7], l[8]) + '\n') + result_pka_file_stripped.write("%3s %3s %s %6s %3s %5s %3s %4s %s" % (l[0], l[1], l[2], l[3], l[4], l[5], l[6], l[7], l[8]) + '\n') + result_pka_file_stripped.close() + result_log.close() + # Now handle the bonds. We have to delete dublicates first. + bonds.sort() + last = bonds[-1] + for i in range(len(bonds) - 2, -1, -1): + if last == bonds[i]: + del bonds[i] + else: + last = bonds[i] + # Now make a selection for known residue + bonds_selected = [] + bonds_ligands = [] + for l in bonds: + if l[0][6:] in ['ASP', 'GLU', 'ARG', 'LYS', 'HIS', 'CYS', 'TYR', 'C-', 'N+']: + bonds_selected.append(l) + else: + bonds_ligands.append(l) + # And now sort it. + bonds_selected.sort(key=lambda residue: int(residue[1])) + # Now write it to file + bonddic = {'=': ' ', ':': ' ', ',': ' ', "'": " "} + for l in bonds_selected: + nb = replace_all(l[7], bonddic) + result_pka_file_bonds.write("%3s %3s %s %7s %7s %9s %17s %s" % (l[0][6:], l[1], l[2][:1], l[3][8:], l[4], l[5], l[6], nb) + '\n') + for l in bonds_ligands: + nb = replace_all(l[7], bonddic) + result_pka_file_bonds.write("%3s %3s %s %7s %7s %9s %17s %s" % (l[0][6:], l[1], l[2][:1], l[3][8:], l[4], l[5], l[6], nb) + '\n') + result_pka_file_bonds.close() + return(result_pka_pkafile) +if runningpymol != 'no': + cmd.extend("getpropka", getpropka) + + +def openpymolfiles(pkafile): + result_pka_pymol_name = pkafile.replace(".pka", ".pml") + result_pka_pymol = open(result_pka_pymol_name, "w") + return(result_pka_pymol, result_pka_pymol_name) + + +def printpropkaresult(list_results, resi, resi_range, resn, resn_range, showresult, ligands_results): + for l in list_results: + if resi != "0" and int(l[1]) in resi_range: + if showresult != 'yes': + print("%3s %3s %s %6s %3s %5s %3s %4s %s" % (l[0], l[1], l[2], l[3], l[4], l[5], l[6], l[7], l[8])) + if resn != "NIL" and l[0] in resn_range and int(l[1]) not in resi_range: + if showresult != 'yes': + print("%3s %3s %s %6s %3s %5s %3s %4s %s" % (l[0], l[1], l[2], l[3], l[4], l[5], l[6], l[7], l[8])) + if showresult == 'yes': + print("%3s %3s %s %6s %3s %5s %3s %4s %s" % (l[0], l[1], l[2], l[3], l[4], l[5], l[6], l[7], l[8])) + for l in ligands_results: + if resn != "NIL" and l[0] in resn_range: + if showresult != 'yes': + print("%3s %3s %s %6s %3s %5s %3s %4s %s" % (l[0], l[1], l[2], l[3], l[4], l[5], l[6], l[7], l[8])) + if showresult == 'yes': + print("%3s %3s %s %6s %3s %5s %3s %4s %s" % (l[0], l[1], l[2], l[3], l[4], l[5], l[6], l[7], l[8])) + + +def importpropkaresult(result_pka_pkafile): + result_pka_file = open(result_pka_pkafile, "r") + list_results = [] + ligands_results = [] + ##bonding_partners = [] + for l in result_pka_file: + if not l.strip(): + continue + else: + # To search for the right lines + if l.strip().split()[0] in ['ASP', 'GLU', 'ARG', 'LYS', 'HIS', 'CYS', 'TYR', 'C-', 'N+'] and len(l.strip().split()) > 20: + list_results.append([l.strip().split()[0], l.strip().split()[1], l.strip().split()[2], l.strip().split()[3], l.strip().split()[4], l.strip().split()[6], l.strip().split()[7], l.strip().split()[8], l.strip().split()[9]]) + # bonding_partners.append(l.strip().split()[11]);bonding_partners.append(l.strip().split()[15]);bonding_partners.append(l.strip().split()[19]) + if l.strip().split()[0] not in ['ASP', 'GLU', 'ARG', 'LYS', 'HIS', 'CYS', 'TYR', 'C-', 'N+'] and len(l.strip().split()) > 20: + ligands_results.append([l.strip().split()[0], l.strip().split()[1], l.strip().split()[2], l.strip().split()[3], l.strip().split()[4], l.strip().split()[6], l.strip().split()[7], l.strip().split()[8], l.strip().split()[9]]) + # bonding_partners.append(l.strip().split()[11]);bonding_partners.append(l.strip().split()[15]);bonding_partners.append(l.strip().split()[19]) + # Sort the result after the residue number and then chain. + list_results.sort(key=lambda residue: int(residue[1])) + list_results.sort(key=lambda chain: chain[2]) + # bonding_partners=uniqifi(bonding_partners) + ##bonding_partners[:] = [x for x in bonding_partners if x != "XXX"] + result_pka_file.close() + return(list_results, ligands_results) + + +def importpropkabonds(result_pka_pkafile): + bonds = [] + result_pka_file_bonds = open(result_pka_pkafile[:-4] + ".bonds", "r") + for l in result_pka_file_bonds: + bonds.append(l.split()) + result_pka_file_bonds.close() + return(bonds) + + +def createdirs(): + if platform.system() == 'Windows': + Newdir = os.getcwd() + "\Results_propka\\" + if platform.system() == 'Linux': + Newdir = os.getcwd() + "/Results_propka/" + if not os.path.exists(Newdir): + os.makedirs(Newdir) + return(Newdir) + + +def openfiles(Newdir, filename, logtime, source): + if source == "upload": + result_pka_pkafile = filename.replace(".pdb", ".pka") + result_pka_input_pkafile = filename.replace(".pdb", ".propka_input") + result_log_name = "%s_Results.log" % (Newdir) + result_pka_file_stripped_name = filename.replace(".pdb", ".stripped") + result_pka_file_bonds_name = filename.replace(".pdb", ".bonds") + if source == "ID": + result_pka_pkafile = "%s%s%s.pka" % (Newdir, filename, logtime) + result_pka_input_pkafile = "%s%s%s.propka_input" % (Newdir, filename, logtime) + result_log_name = "%s_Results.log" % (Newdir) + result_pka_file_stripped_name = "%s%s%s.stripped" % (Newdir, filename, logtime) + result_pka_file_bonds_name = "%s%s%s.bonds" % (Newdir, filename, logtime) + if platform.system() == 'Windows': + filepath = "\\" + if platform.system() == 'Linux': + filepath = "/" + # Open the files + result_pka_file = open(result_pka_pkafile, "w") + result_input_pka_file = open(result_pka_input_pkafile, "w") + result_log = open(result_log_name, "a") + result_pka_file_stripped = open(result_pka_file_stripped_name, "w") + result_pka_file_bonds = open(result_pka_file_bonds_name, "w") + return(result_pka_file, result_input_pka_file, result_log, filepath, result_pka_pkafile, result_pka_file_stripped, result_pka_file_bonds) + + +def ResiRange(resi): + resi = resi.split('.') + resiList = [] + for i in resi: + if '-' in i: + tmp = i.split('-') + resiList.extend(list(range(int(tmp[0]), int(tmp[-1]) + 1))) + if '-' not in i: + resiList.append(int(i)) + return(resiList) + + +def ResnRange(resn): + resn_split = resn.split('.') + resn_range = [resnr.upper() for resnr in resn_split] + return(resn_range) + + +def ChainRange(chain): + chainstring = chain.replace(".", "+").upper() + return(chainstring) + + +def writepymolcmd(newmolecule, pkafile, verbose, makebonds): + list_results, ligands_results = importpropkaresult(pkafile) + # Now find the available bonding partners that pymol knows of + bonding_partners = [] + bonding_partners_str = cmd.get_pdbstr("%s and resn * and not resn ASP+GLU+ARG+LYS+HIS+CYS+TYR+GLN+ASN+SER+THR+GLY+PHE+LEU+ALA+ILE+TRP+MET+PRO+VAL+HOH" % (newmolecule)) + for i in range(len(bonding_partners_str.splitlines()) - 1): + bonding_partners_split = bonding_partners_str.splitlines()[i].split() + if bonding_partners_split[0] == "HETATM" or bonding_partners_split[0] == "ATOM": + bonding_partners_single = bonding_partners_split[3] + bonding_partners.append(bonding_partners_single) + bonding_partners = uniqifi(bonding_partners) + if verbose == 'yes': + print("And other possible bonding partners is: %s" % bonding_partners) + # Read in the bond file, if it exists + writebonds = "no" + if os.path.isfile(pkafile[:-4] + ".bonds") and makebonds == "yes": + bonds = importpropkabonds(pkafile) + writebonds = "yes" + # Open the pymol command file for writing + files_pka_pymol = openpymolfiles(pkafile) + result_pka_pymol = files_pka_pymol[0] + result_pka_pymol_name = files_pka_pymol[1] + # Make some dictionary for propka->pymol name conversion + dictio = {'ASP': 'CG', 'GLU': 'CD', 'ARG': 'CZ', 'LYS': 'NZ', 'HIS': 'CG', 'CYS': 'SG', 'TYR': 'OH', 'C-': 'C', 'N+': 'N', 'NTR': 'N', 'CTR': 'C', 'GLN': 'CD', 'ASN': 'CG', 'SER': 'OG', 'THR': 'OG1', 'GLY': 'CA', 'PHE': 'CZ', 'LEU': 'CG', 'ALA': 'CB', 'ILE': 'CD1', 'TRP': 'NE1', 'MET': 'SD', 'PRO': 'CG', 'VAL': 'CB'} + dictio2 = {'ASP': 'D', 'GLU': 'E', 'ARG': 'R', 'LYS': 'K', 'HIS': 'H', 'CYS': 'C', 'TYR': 'Y', 'C-': 'C-', 'N+': 'N+'} + # This list is from: http://en.wikipedia.org/wiki/Protein_pKa_calculations + pkaaminoacid = ['ASP', 'GLU', 'ARG', 'LYS', 'HIS', 'CYS', 'TYR'] + pkadictio = {'ASP': 3.9, 'GLU': 4.3, 'ARG': 12.0, 'LYS': 10.5, 'HIS': 6.0, 'CYS': 8.3, 'TYR': 10.1} + # Now start write to the file. + # Try to make silent + # result_pka_pymol.write("cmd.feedback('disable','all','actions')\n") + # result_pka_pymol.write("cmd.feedback('disable','all','results')\n") + # Change the GUI width, to make the long names possible. + result_pka_pymol.write("cmd.set('internal_gui_width','360')\n") + # Set fonts + result_pka_pymol.write("cmd.set('label_font_id','12')\n") + result_pka_pymol.write("cmd.set('label_size','-0.5')\n") + result_pka_pymol.write("cmd.set('label_color','grey')\n") + # No auto zoom the new objects + result_pka_pymol.write("cmd.set('auto_zoom','off')\n") + # The name for the molecules are defined here + pkamolecule = "%spKa" % (newmolecule) + pkalabelmolecule = "%sLab" % (newmolecule) + # Create the groups now, so they come in order. They will be empty + result_pka_pymol.write("cmd.group('%sResi','Res*')\n" % (newmolecule)) + result_pka_pymol.write("cmd.group('%sLigands','Lig*')\n" % (newmolecule)) + if writebonds == "yes": + result_pka_pymol.write("cmd.group('%sBonds','%sBond*')\n" % (newmolecule, newmolecule)) + # Create new empty pymol pka molecules. For pka atoms and its label. This is a "bucket" we where we will put in the atoms together. + result_pka_pymol.write("cmd.create('%s','None')\n" % (pkamolecule)) + result_pka_pymol.write("cmd.create('%s','None')\n" % (pkalabelmolecule)) + # Now make the pka atoms and alter, color and such + for l in list_results: + name = dictio[l[0]] + resn = dictio2[l[0]] + resi = l[1] + chain = l[2] + pka = l[3] + buried = l[4] + if "*" in pka: + pka = pka.replace("*", "") + comment = "*Coupled residue" + else: + comment = "" + if l[0] in pkaaminoacid: + pkadiff = (float(pka) - pkadictio[l[0]]) + pkadiff = "(%s)" % pkadiff + if pka == "99.99": + pkadiff = "" + else: + pkadiff = "" + # Make the selection for which atom to copy + newselection = ("/%s//%s/%s/%s" % (newmolecule, chain, resi, name)) + protselect = ("%sRes_%s%s%s" % (newmolecule, chain, resn, resi)) + result_pka_pymol.write("cmd.select('%s','byres %s')\n" % (protselect, newselection)) + result_pka_pymol.write("cmd.show('sticks','byres %s')\n" % (protselect)) + # The temporary name + tempname = ("%s%s%s%s" % (pkamolecule, chain, resi, name)) + tempnamelabel = ("%s%s%s%s" % (pkalabelmolecule, chain, resi, name)) + tempselect = ("/%s//%s/%s" % (tempname, chain, resi)) + tempselectlabel = ("/%s//%s/%s" % (tempnamelabel, chain, resi)) + # Copy the atom, call it by the residue name + result_pka_pymol.write("cmd.create('%s','%s',quiet=1)\n" % (tempname, newselection)) + # Alter the name and the b value of the newly created atom + result_pka_pymol.write("cmd.alter('%s','b=%s')\n" % (tempselect, pka)) + result_pka_pymol.write("cmd.alter('%s','vdw=0.5')\n" % (tempselect)) + result_pka_pymol.write("cmd.alter('%s','name=%s%s%s')\n" % (tempselect, '"', pka, '"')) + # Now create a fake label atom, and translate it + result_pka_pymol.write("cmd.create('%s','%s',quiet=1)\n" % (tempnamelabel, tempname)) + movelabelxyz = (1.5, 0, 0) + result_pka_pymol.write("cmd.translate('[%s,%s,%s]','%s',camera=0)\n" % (movelabelxyz[0], movelabelxyz[1], movelabelxyz[2], tempnamelabel)) + # Labelling alternate positions are not allowed, so we delete that attribute for the label atoms. + result_pka_pymol.write("cmd.alter('%s','alt=%s%s')\n" % (tempselectlabel, '"', '"')) + result_pka_pymol.write("cmd.label('%s','text_type=%spka=%s%s Bu:%s%s%s%s')\n" % (tempselectlabel, '"', pka, pkadiff, buried, '%', comment, '"')) + # Now put the atoms into a bucket of atoms + result_pka_pymol.write("cmd.create('%s','%s or (%s)',quiet=1)\n" % (pkamolecule, pkamolecule, tempselect)) + result_pka_pymol.write("cmd.create('%s','%s or (%s)',quiet=1)\n" % (pkalabelmolecule, pkalabelmolecule, tempselectlabel)) + # Remove the temporary atoms + result_pka_pymol.write("cmd.remove('%s')\n" % (tempname)) + result_pka_pymol.write("cmd.remove('%s')\n" % (tempnamelabel)) + # Delete the temporary molecule/selection + result_pka_pymol.write("cmd.delete('%s')\n" % (tempname)) + result_pka_pymol.write("cmd.delete('%s')\n" % (tempnamelabel)) + # Group the resi together + result_pka_pymol.write("cmd.group('%sResi','%sRes*')\n" % (newmolecule, newmolecule)) + for l in ligands_results: + resn = l[0] + atom = l[1] + chain = l[2] + pka = l[3] + buried = l[4] + if verbose == 'yes': + print("Ligand. resn:%s atom:%s chain:%s pka:%s buried:%s" % (resn, atom, chain, pka, buried)) + if Check_bonding_partners(bonding_partners, resn)[0]: + if "*" in pka: + pka = pka.replace("*", "") + comment = "*Coupled residue" + else: + comment = "" + # Make the selection for which atom to copy + ligselection = ("/%s and chain %s and resn %s and name %s" % (newmolecule, chain, resn, atom)) + ligselect = ("%sLig_%s%s%s" % (newmolecule, chain, resn, atom)) + result_pka_pymol.write("cmd.select('%s','%s')\n" % (ligselect, ligselection)) + result_pka_pymol.write("cmd.show('sticks','byres %s')\n" % (ligselect)) + result_pka_pymol.write("cmd.util.cbap('byres %s')\n" % (ligselect)) + # The temporary name + tempname = ("%s%s%s%s" % (pkamolecule, chain, resn, atom)) + tempnamelabel = ("%s%s%s%s" % (pkalabelmolecule, chain, resn, atom)) + tempselect = ("/%s and chain %s and resn %s" % (tempname, chain, resn)) + tempselectlabel = ("/%s and chain %s and resn %s" % (tempnamelabel, chain, resn)) + # Copy the atom, call it by the residue name + result_pka_pymol.write("cmd.create('%s','%s',quiet=1)\n" % (tempname, ligselection)) + # Alter the name and the b value of the newly created atom + result_pka_pymol.write("cmd.alter('%s','b=%s')\n" % (tempselect, pka)) + result_pka_pymol.write("cmd.alter('%s','vdw=0.5')\n" % (tempselect)) + result_pka_pymol.write("cmd.alter('%s','name=%s%s%s')\n" % (tempselect, '"', pka, '"')) + # Now create a fake label atom, and translate it + result_pka_pymol.write("cmd.create('%s','%s',quiet=1)\n" % (tempnamelabel, tempname)) + movelabelxyz = (1.5, 0, 0) + result_pka_pymol.write("cmd.translate('[%s,%s,%s]','%s',camera=0)\n" % (movelabelxyz[0], movelabelxyz[1], movelabelxyz[2], tempnamelabel)) + # Labelling alternate positions are not allowed, so we delete that attribute for the label atoms. + result_pka_pymol.write("cmd.alter('%s','alt=%s%s')\n" % (tempselectlabel, '"', '"')) + result_pka_pymol.write("cmd.label('%s','text_type=%spka=%s Bu:%s%s%s%s')\n" % (tempselectlabel, '"', pka, buried, '%', comment, '"')) + # Now put the atoms into a bucket of atoms + result_pka_pymol.write("cmd.create('%s','%s or (%s)',quiet=1)\n" % (pkamolecule, pkamolecule, tempselect)) + result_pka_pymol.write("cmd.create('%s','%s or (%s)',quiet=1)\n" % (pkalabelmolecule, pkalabelmolecule, tempselectlabel)) + # Remove the temporary atoms + result_pka_pymol.write("cmd.remove('%s')\n" % (tempname)) + result_pka_pymol.write("cmd.remove('%s')\n" % (tempnamelabel)) + # Delete the temporary molecule/selection + result_pka_pymol.write("cmd.delete('%s')\n" % (tempname)) + result_pka_pymol.write("cmd.delete('%s')\n" % (tempnamelabel)) + # Group the resi together + result_pka_pymol.write("cmd.group('%sLigands','%sLig*')\n" % (newmolecule, newmolecule)) + # Finish the pka atoms, and show spheres + result_pka_pymol.write("cmd.show('spheres','%s')\n" % (pkamolecule)) + result_pka_pymol.write("cmd.spectrum('b','red_white_blue',selection='%s',minimum='0',maximum='14')\n" % (pkamolecule)) + result_pka_pymol.write("cmd.alter('%s and name 99.9','vdw=0.8')\n" % (pkamolecule)) + result_pka_pymol.write("cmd.show('spheres','%s and name 99.9')\n" % (pkamolecule)) + result_pka_pymol.write("cmd.color('sulfur','%s and name 99.9')\n" % (pkamolecule)) + # Now we make the bonds + if writebonds == "yes": + Bondgroups = [] + naturalaminoacids = ['ASP', 'GLU', 'ARG', 'LYS', 'HIS', 'CYS', 'TYR', 'NTR', 'N+', 'CTR', 'C-', 'GLN', 'ASN', 'SER', 'THR', 'GLY', 'PHE', 'LEU', 'ALA', 'ILE', 'TRP', 'MET', 'PRO', 'VAL'] + for l in bonds: + if l[0] in naturalaminoacids: + name = dictio[l[0]] + resi = l[1] + chain = l[2] + desolvation = l[6][12:] + pkachange = l[11] + NBresi = l[8][3:] + NBchain = l[9] + NBbond = l[-1][:2] + if l[8][:3] in naturalaminoacids: + NBname, cutoff = BondTypeName(dictio[l[8][:3]], NBbond) + fromselection = ("/%s//%s/%s/%s" % (newmolecule, chain, resi, name)) + toselection = ("/%s//%s/%s/%s" % (newmolecule, NBchain, NBresi, NBname)) + if l[8][:3] == 'NTR': + extind = cmd.identify("chain %s and name N" % (NBchain))[0] + toselection = ("/%s and id %s and name N" % (newmolecule, extind)) + NBresi = "N+" + if l[8][:3] == 'CTR': + extind = cmd.identify("chain %s and name C" % (NBchain))[-1] + toselection = ("/%s and id %s and name C" % (newmolecule, extind)) + NBresi = "C-" + distname = ("%s_%s%s%s%s_%s_%s" % (newmolecule, chain, resi, NBchain, NBresi, NBbond, pkachange)) + result_pka_pymol.write("cmd.distance('%s','%s','%s'%s)\n" % (distname, fromselection, toselection, cutoff)) + result_pka_pymol.write("cmd.color('%s','%s')\n" % (SetDashColor(NBbond), distname)) + # result_pka_pymol.write("cmd.disable('%s')\n"%(distname)) + Bondgroups.append("%s%s" % (chain, resi)) + if l[8][:3] not in naturalaminoacids and Check_bonding_partners(bonding_partners, l[8])[0]: + cutoff = "" + NBresn = Check_bonding_partners(bonding_partners, l[8])[1] + NBname = l[8][len(NBresn):] + "*" + fromselection = ("/%s//%s/%s/%s" % (newmolecule, chain, resi, name)) + toselection = ("/%s and chain %s and resn %s and name %s" % (newmolecule, NBchain, NBresn, NBname)) + if verbose == 'yes': + print("Res->Ligand: (%s) -> (%s)" % (fromselection, toselection)) + result_pka_pymol.write("cmd.show('sticks','%s')\n" % (toselection)) + distname = ("%s_%s%s%s_%s_%s" % (newmolecule, chain, resi, NBresn, NBbond, pkachange)) + result_pka_pymol.write("cmd.distance('%s','%s','%s'%s)\n" % (distname, fromselection, toselection, cutoff)) + result_pka_pymol.write("cmd.color('%s','%s')\n" % (SetDashColor(NBbond), distname)) + # result_pka_pymol.write("cmd.disable('%s')\n"%(distname)) + Bondgroups.append("%s%s" % (chain, resi)) + if l[0] in bonding_partners: + resn = l[0] + atom = l[1] + chain = l[2] + desolvation = l[6][12:] + pkachange = l[11] + NBresi = l[8][3:] + NBchain = l[9] + NBbond = l[-1][:2] + if not Check_bonding_partners(bonding_partners, l[8])[0]: + NBname, cutoff = BondTypeName(dictio[l[8][:3]], NBbond) + fromselection = ("/%s and chain %s and resn %s and name %s" % (newmolecule, chain, resn, atom)) + toselection = ("/%s//%s/%s/%s" % (newmolecule, NBchain, NBresi, NBname)) + if l[8][:3] == 'NTR': + extind = cmd.identify("chain %s and name N" % (NBchain))[0] + toselection = ("/%s and id %s and name N" % (newmolecule, extind)) + NBresi = "N+" + if l[8][:3] == 'CTR': + extind = cmd.identify("chain %s and name C" % (NBchain))[-1] + toselection = ("/%s and id %s and name C" % (newmolecule, extind)) + NBresi = "C-" + distname = ("%s_%s%s%s%s%s_%s_%s" % (newmolecule, chain, resn, atom, NBchain, NBresi, NBbond, pkachange)) + result_pka_pymol.write("cmd.distance('%s','%s','%s'%s)\n" % (distname, fromselection, toselection, cutoff)) + result_pka_pymol.write("cmd.color('%s','%s')\n" % (SetDashColor(NBbond), distname)) + Bondgroups.append("%s%s%s" % (chain, resn, atom)) + # result_pka_pymol.write("cmd.disable('%s')\n"%(distname)) + if Check_bonding_partners(bonding_partners, l[8])[0]: + cutoff = "" + NBresn = Check_bonding_partners(bonding_partners, l[8])[1] + NBname = l[8][len(NBresn):] + "*" + fromselection = ("/%s and chain %s and resn %s and name %s" % (newmolecule, chain, resn, atom)) + toselection = ("/%s and chain %s and resn %s and name %s" % (newmolecule, NBchain, NBresn, NBname)) + if verbose == 'yes': + print("Ligand->Ligand: (%s) -> (%s)" % (fromselection, toselection)) + result_pka_pymol.write("cmd.show('sticks','%s')\n" % (toselection)) + distname = ("%s_%s%s%s%s_%s_%s" % (newmolecule, chain, resn, atom, NBresn, NBbond, pkachange)) + result_pka_pymol.write("cmd.distance('%s','%s','%s'%s)\n" % (distname, fromselection, toselection, cutoff)) + result_pka_pymol.write("cmd.color('%s','%s')\n" % (SetDashColor(NBbond), distname)) + # result_pka_pymol.write("cmd.disable('%s')\n"%(distname)) + Bondgroups.append("%s%s%s" % (chain, resn, atom)) + Bondgroups = uniqifi(Bondgroups) + for l in Bondgroups: + result_pka_pymol.write("cmd.group('%sBonds_%s','%s_%s*')\n" % (newmolecule, l, newmolecule, l)) + result_pka_pymol.write("cmd.disable('%sBonds_%s')\n" % (newmolecule, l)) + result_pka_pymol.write("cmd.group('%sBonds','%sBonds_*')\n" % (newmolecule, newmolecule)) + result_pka_pymol.write("cmd.set('auto_zoom','on')\n") + # result_pka_pymol.write("cmd.feedback('enable','all','actions')\n") + # result_pka_pymol.write("cmd.feedback('enable','all','results')\n") + result_pka_pymol.close() + return(result_pka_pymol_name) + + +def replace_all(text, dic): + for i, j in dic.items(): + text = text.replace(i, j) + return(text) + + +def uniqifi(seq, idfun=None): + # Order preserving + if idfun is None: + def idfun(x): + return x + seen = {} + result = [] + for item in seq: + marker = idfun(item) + if marker in seen: + continue + seen[marker] = 1 + result.append(item) + return(result) + + +def BondTypeName(NBname, NBbond): + if NBbond == "SH": + cutoff = "" + return(NBname, cutoff) + if NBbond == "BH": + cutoff = "" + return("N", cutoff) + else: + cutoff = "" + return(NBname, cutoff) + + +def Check_bonding_partners(bonding_partners, NBname): + answer = False + for l in bonding_partners: + if l in NBname: + answer = True + NBname = l + break + else: + answer = False + return(answer, NBname) + + +def SetDashColor(NBbond): + if NBbond == "SH": + color = "brightorange" + if NBbond == "BH": + color = "lightorange" + if NBbond == "CC": + color = "red" + return(color) diff --git a/.pymol/startup/propka.pyc b/.pymol/startup/propka.pyc new file mode 100644 index 0000000..6ee65c3 Binary files /dev/null and b/.pymol/startup/propka.pyc differ diff --git a/.pymol/startup/show_contact_map.py b/.pymol/startup/show_contact_map.py new file mode 100644 index 0000000..ffa2be1 --- /dev/null +++ b/.pymol/startup/show_contact_map.py @@ -0,0 +1,58 @@ +#!/usr/bin/env python + +import tkFileDialog +import tkSimpleDialog +from pymol import cmd +import math +from distutils.util import strtobool + +def __init__(self): + self.menuBar.addmenuitem('Plugin', 'command', 'Contact map', label='Show contact map', command=lambda s=self: load_map(s)) + +def load_map(app): + contacts_file = tkFileDialog.askopenfile(initialdir=".", title="Load contact map") + if contacts_file is None: + return + show_contact_map(file=contacts_file) + +def show_contact_map(selection="all", file=None, as_bonds=1, quiet=0): + """ Visualize contact map + """ + as_bonds = bool(as_bonds) + quiet = bool(quiet) + contacts = read_contacts(file) + + for contact in contacts: + name = "contact{}-{}".format(contact[0], contact[1]) + seleA = "id %s and %s" % (contact[0], selection) + seleB = "id %s and %s" % (contact[1], selection) + + # skip non-exsistant + if cmd.select("seleA", seleA) == "0" or cmd.select("seleB", seleB) == "0": + print("%s to %s not found!" % (seleA, seleB)) + + if not quiet: + print("New bond: residue %s to %s" % (contact[0], contact[1])) + + if as_bonds: + cmd.bond(seleA, seleB) + else: + cmd.distance(name, seleA, seleB) + + cmd.delete('seleA') + cmd.delete('seleB') + + +def read_contacts(contacts_file, separator=' '): + contacts = [] + if not isinstance(contacts_file, file): + contacts_file = open(contacts_file, 'r') + for line in contacts_file.readlines(): + split = line.split(",") + resA = int(split[0]) + resB = int(split[1]) + contacts.append([resA, resB]) + contacts_file.close() + return contacts + +cmd.extend("contact_map", show_contact_map) diff --git a/.pymol/startup/show_contact_map.pyc b/.pymol/startup/show_contact_map.pyc new file mode 100644 index 0000000..8574cc6 Binary files /dev/null and b/.pymol/startup/show_contact_map.pyc differ diff --git a/.pymol/startup/spectrum_states.py b/.pymol/startup/spectrum_states.py new file mode 100644 index 0000000..ade94b3 --- /dev/null +++ b/.pymol/startup/spectrum_states.py @@ -0,0 +1,92 @@ +''' +http://pymolwiki.org/index.php/spectrum_states + +(c) 2011 Takanori Nakane and Thomas Holder + +License: BSD-2-Clause +''' + +from __future__ import print_function + +from pymol import cmd, CmdException + + +def spectrum_states(selection='all', representations='cartoon ribbon', + color_list='blue cyan green yellow orange red', + first=1, last=0, quiet=1): + ''' +DESCRIPTION + + Color each state in a multi-state object different. + +USAGE + + spectrum_states [ selection [, representations [, color_list [, first [, last ]]]]] + +ARGUMENTS + + selection = string: object names (works with complete objects only) + {default: all} + + representations = string: space separated list of representations + {default: cartoon ribbon} + + color_list = string: space separated list of colors {default: blue cyan + green yellow orange red} + +SEE ALSO + + spectrum, spectrumany + ''' + from math import floor, ceil + + first, last, quiet = int(first), int(last), int(quiet) + colors = color_list.split() + if len(colors) < 2: + print(' Error: please provide at least 2 colors') + raise CmdException + + colvec = [cmd.get_color_tuple(i) for i in colors] + + # filter for valid _color settings + settings = [] + for r in representations.split(): + if r[-1] == 's': + r = r[:-1] + s = r + '_color' + if s in cmd.setting.name_list: + settings.append(s) + elif not quiet: + print(' Warning: no such setting:', s) + + # object names only + selection = ' '.join(cmd.get_object_list('(' + selection + ')')) + if cmd.count_atoms(selection) == 0: + print(' Error: empty selection') + raise CmdException + + if last < 1: + last = cmd.count_states(selection) + + val_range = int(last - first + 1) + if val_range < 2: + print(' Error: no spectrum possible, need more than 1 state') + raise CmdException + + for i in range(val_range): + p = float(i) / (val_range - 1) * (len(colvec) - 1) + p0, p1 = int(floor(p)), int(ceil(p)) + ii = (p - p0) + col_list = [colvec[p1][j] * ii + colvec[p0][j] * (1.0 - ii) for j in range(3)] + col_name = '0x%02x%02x%02x' % (col_list[0] * 255, col_list[1] * 255, col_list[2] * 255) + for s in settings: + cmd.set(s, col_name, selection, state=i + first) + +cmd.extend('spectrum_states', spectrum_states) + +# tab-completion of arguments +cmd.auto_arg[0]['spectrum_states'] = cmd.auto_arg[0]['disable'] +cmd.auto_arg[1]['spectrum_states'] = [cmd.auto_arg[0]['show'][0], 'representation', ' '] +cmd.auto_arg[2]['spectrum_states'] = [cmd.auto_arg[0]['color'][0], 'color', ' '] + +# vi:expandtab:smarttab diff --git a/.pymol/startup/spectrum_states.pyc b/.pymol/startup/spectrum_states.pyc new file mode 100644 index 0000000..4491a2c Binary files /dev/null and b/.pymol/startup/spectrum_states.pyc differ diff --git a/.pymol/startup/tmalign.py b/.pymol/startup/tmalign.py new file mode 100644 index 0000000..46d5257 --- /dev/null +++ b/.pymol/startup/tmalign.py @@ -0,0 +1,257 @@ +''' +(c) 2011-2012 Thomas Holder, MPI for Developmental Biology +''' + +from __future__ import print_function + +__author__ = 'Thomas Holder' +__version__ = '1.1' +__license__ = 'BSD-2-Clause' + +from pymol import cmd, CmdException + + +def save_pdb_without_ter(filename, selection, **kwargs): + ''' +DESCRIPTION + +Save PDB file without TER records. External applications like TMalign and +DynDom stop reading PDB files at TER records, which might be undesired in +case of missing loops. + ''' + v = cmd.get_setting_boolean('pdb_use_ter_records') + if v: + cmd.unset('pdb_use_ter_records') + cmd.save(filename, selection, **kwargs) + if v: + cmd.set('pdb_use_ter_records') + + +def alignwithanymethod(mobile, target, methods='align super cealign tmalign', + async=1, quiet=1): + ''' +DESCRIPTION + +Align copies of mobile to target with several alignment methods + +ARGUMENTS + +mobile = string: atom selection + +target = string: atom selection + +methods = string: space separated list of PyMOL commands which take +arguments "mobile" and "target" (in any order) {default: align super +cealign tmalign} + ''' + import threading + import time + methods = methods.split() + async, quiet = int(async), int(quiet) + mobile_obj = cmd.get_object_list('first (' + mobile + ')')[0] + + def myalign(method): + newmobile = cmd.get_unused_name(mobile_obj + '_' + method) + cmd.create(newmobile, mobile_obj) + start = time.time() + cmd.do('%s mobile=%s in %s, target=%s' % (method, newmobile, mobile, target)) + if not quiet: + print('Finished: %s (%.2f sec)' % (method, time.time() - start)) + + for method in methods: + if async: + t = threading.Thread(target=myalign, args=(method,)) + t.setDaemon(1) + t.start() + else: + myalign(method) + + +def tmalign(mobile, target, args='', exe='TMalign', ter=0, transform=1, object=None, quiet=0): + ''' +DESCRIPTION + +TMalign wrapper + +Reference: Y. Zhang and J. Skolnick, Nucl. Acids Res. 2005 33, 2302-9 +http://zhanglab.ccmb.med.umich.edu/TM-align/ + +USAGE + +tmalign mobile, target [, args [, exe ]] + +ARGUMENTS + +mobile, target = string: atom selections + +args = string: Extra arguments like -d0 5 -L 100 + +exe = string: Path to TMalign executable {default: TMalign} + +ter = 0/1: If ter=0, then ignore chain breaks because TMalign will stop +at first TER record {default: 0} + +SEE ALSO + +tmscore, mmalign + ''' + import subprocess + import tempfile + import os + import re + + ter, quiet = int(ter), int(quiet) + + mobile_filename = tempfile.mktemp('.pdb', 'mobile') + target_filename = tempfile.mktemp('.pdb', 'target') + matrix_filename = tempfile.mktemp('.txt', 'matrix') + mobile_ca_sele = '(%s) and (not hetatm) and name CA and alt +A' % (mobile) + target_ca_sele = '(%s) and (not hetatm) and name CA and alt +A' % (target) + + if ter: + save = cmd.save + else: + save = save_pdb_without_ter + save(mobile_filename, mobile_ca_sele) + save(target_filename, target_ca_sele) + + exe = cmd.exp_path(exe) + args = [exe, mobile_filename, target_filename, '-m', matrix_filename] + args.split() + + try: + process = subprocess.Popen(args, stdout=subprocess.PIPE) + lines = process.stdout.readlines() + except OSError: + print('Cannot execute "%s", please provide full path to TMscore or TMalign executable' % (exe)) + raise CmdException + finally: + os.remove(mobile_filename) + os.remove(target_filename) + + # TMalign >= 2012/04/17 + if os.path.exists(matrix_filename): + lines += open(matrix_filename).readlines() + os.remove(matrix_filename) + + r = None + re_score = re.compile(r'TM-score\s*=\s*(\d*\.\d*)') + rowcount = 0 + matrix = [] + line_it = iter(lines) + alignment = [] + for line in line_it: + if 4 >= rowcount > 0: + if rowcount >= 2: + a = list(map(float, line.split())) + matrix.extend(a[2:5]) + matrix.append(a[1]) + rowcount += 1 + elif line.lower().startswith(' -------- rotation matrix'): + rowcount = 1 + elif line.startswith('(":" denotes'): + alignment = [line_it.next().rstrip() for i in range(3)] + else: + match = re_score.search(line) + if match is not None: + r = float(match.group(1)) + if not quiet: + print(line.rstrip()) + + if not quiet: + for i in range(0, len(alignment[0]) - 1, 78): + for line in alignment: + print(line[i:i + 78]) + print('') + + assert len(matrix) == 3 * 4 + matrix.extend([0, 0, 0, 1]) + + if int(transform): + cmd.transform_selection('byobject (%s)' % (mobile), matrix, homogenous=1) + + # alignment object + if object is not None: + mobile_idx, target_idx = [], [] + space = {'mobile_idx': mobile_idx, 'target_idx': target_idx} + cmd.iterate(mobile_ca_sele, 'mobile_idx.append("%s`%d" % (model, index))', space=space) + cmd.iterate(target_ca_sele, 'target_idx.append("%s`%d" % (model, index))', space=space) + for i, aa in enumerate(alignment[0]): + if aa == '-': + mobile_idx.insert(i, None) + for i, aa in enumerate(alignment[2]): + if aa == '-': + target_idx.insert(i, None) + if (len(mobile_idx) == len(target_idx) == len(alignment[2])): + cmd.rms_cur( + ' '.join(idx for (idx, m) in zip(mobile_idx, alignment[1]) if m in ':.'), + ' '.join(idx for (idx, m) in zip(target_idx, alignment[1]) if m in ':.'), + cycles=0, matchmaker=4, object=object) + else: + print('Could not load alignment object') + + if not quiet and r is not None: + print('Found in output TM-score = %.4f' % (r)) + + return r + + +def tmscore(mobile, target, args='', exe='TMscore', quiet=0, **kwargs): + ''' +DESCRIPTION + +TMscore wrapper + +Reference: Yang Zhang and Jeffrey Skolnick, Proteins 2004 57: 702-710 +http://zhanglab.ccmb.med.umich.edu/TM-score/ + +ARGUMENTS + +mobile, target = string: atom selections + +args = string: Extra arguments like -d 5 + +exe = string: Path to TMscore executable {default: TMscore} + +ter = 0/1: If ter=0, then ignore chain breaks because TMscore will stop +at first TER record {default: 0} + +SEE ALSO + +tmalign, mmalign + ''' + kwargs.pop('_self', None) + return tmalign(mobile, target, args, exe, quiet=quiet, **kwargs) + + +def mmalign(mobile, target, args='', exe='MMalign', ter=0, transform=1, quiet=0): + ''' +DESCRIPTION + +MMalign wrapper + +Reference: S. Mukherjee and Y. Zhang, Nucleic Acids Research 2009; 37: e83 +http://zhanglab.ccmb.med.umich.edu/MM-align/ + +SEE ALSO + +tmalign, tmscore + ''' + return tmalign(mobile, target, args, exe, ter, transform, quiet=quiet) + +# pymol commands +cmd.extend('alignwithanymethod', alignwithanymethod) +cmd.extend('tmalign', tmalign) +cmd.extend('tmscore', tmscore) +cmd.extend('mmalign', tmalign) + +# autocompletion +cmd.auto_arg[0].update({ + 'tmalign': cmd.auto_arg[0]['align'], + 'tmscore': cmd.auto_arg[0]['align'], + 'mmalign': cmd.auto_arg[0]['align'], +}) +cmd.auto_arg[1].update({ + 'tmalign': cmd.auto_arg[1]['align'], + 'tmscore': cmd.auto_arg[1]['align'], + 'mmalign': cmd.auto_arg[1]['align'], +}) diff --git a/.pymol/startup/tmalign.pyc b/.pymol/startup/tmalign.pyc new file mode 100644 index 0000000..d2cc90c Binary files /dev/null and b/.pymol/startup/tmalign.pyc differ diff --git a/.pymol/startup/visualize_dca_scores.py b/.pymol/startup/visualize_dca_scores.py new file mode 100644 index 0000000..1862251 --- /dev/null +++ b/.pymol/startup/visualize_dca_scores.py @@ -0,0 +1,87 @@ +#!/usr/bin/env python + +import tkFileDialog +import tkSimpleDialog +from pymol import cmd +import math + +DEFAULT_CUTOFF = 10 +DEFAULT_FLOOR = 0.6 +DEFAULT_MAX_CONTACTS = -1 + +def __init__(self): + self.menuBar.addmenuitem('Plugin', 'command', 'MSA Scores', label='MSA Scores', command=lambda s=self: load_scores_dialog(s)) + +def load_scores_dialog(app): + scores_file = tkFileDialog.askopenfile(initialdir=".", title="Load scores") + if scores_file is None: + return + cutoff = tkSimpleDialog.askinteger("Cutoff", "Backbone cutoff", initialvalue=DEFAULT_CUTOFF, minvalue=1) + if cutoff is None: + return + floor = tkSimpleDialog.askfloat("Minimal Score", "Minimal Score", initialvalue=DEFAULT_FLOOR) + if floor is None: + return + max_contacts = tkSimpleDialog.askfloat("Limit contacts", "Maximal contacts to show", initialvalue=DEFAULT_MAX_CONTACTS) + if max_contacts is None: + return + show_scores(scores=scores_file, cutoff=cutoff, floor=floor, max_contacts=max_contacts) + +def show_scores(selection='all', scores="scores.csv", floor=DEFAULT_FLOOR, cutoff=DEFAULT_CUTOFF, max_contacts=None): + """ Visualize score csv + selection: Selection for visualization + floor: minimal score for visualization + cutoff: minimal distance between amino acids + scores: scores file for visualization + """ + cutoff = int(cutoff) + floor = float(floor) + scores = read_scores(scores) + + max_contacts = int(max_contacts) + if max_contacts is None or max_contacts == -1: + max_contacts = len(scores) + + # filter for cutoff + scores = list(filter(lambda x: x[1] - x[0] >= cutoff, scores)) + + # filter for score floor + scores = list(filter(lambda x: x[2] >= floor, scores)) + + # show distances + counter = 0 + for score in scores: + if counter > max_contacts: + print("max contacts reached! (%s contacts not shown)" % (len(scores)-int(max_contacts))) + break + print("new msa contact between %s and %s" % (score[0], score[1])) + name = "msa{}-{}_{:.3f}".format(score[0], score[1], score[2]) + seleA = "resid %s and n. CA and %s" % (score[0], selection) + seleB = "resid %s and n. CA and %s" % (score[1], selection) + if cmd.select("seleA", seleA) == "0" or cmd.select("seleB", seleB) == "0": + continue + cmd.distance(name, seleA, seleB) + counter += 1 + + cmd.delete('seleA') + cmd.delete('seleB') + + +def read_scores(scores_file): + scores = [] + if not isinstance(scores_file, file): + scores_file = open(scores_file, 'r') + for line in scores_file.readlines(): + split = line.split(",") + resA = int(split[0]) + resB = int(split[1]) + score = float(split[2]) + dist = abs(resB - resA) + scores.append([resA, resB, score, dist]) + scores_file.close() + return scores + +def remove_scores(): + cmd.delete("msa*") +cmd.extend("remove_scores", remove_scores) +cmd.extend("show_scores", show_scores) diff --git a/.pymol/startup/visualize_dca_scores.pyc b/.pymol/startup/visualize_dca_scores.pyc new file mode 100644 index 0000000..1ac3211 Binary files /dev/null and b/.pymol/startup/visualize_dca_scores.pyc differ diff --git a/.pymolrc b/.pymolrc new file mode 100644 index 0000000..4832bb2 --- /dev/null +++ b/.pymolrc @@ -0,0 +1,4 @@ +set opaque_background, 0 +log_open /tmp/pymol_log.pml +set cartoon_side_chain_helper, 1 +set cartoon_highlight_color, grey50