update from laptop

.pymol/startup/anglebetweenhelices.py

@@ -0,0 +1,316 @@
+'''
+(c) 2010 Thomas Holder
+'''
+
+from pymol import cmd, stored, CmdException
+from chempy import cpv
+import math
+
+if cmd.get_version()[1] < 1.2:
+    def get_unused_name(name):
+        import random
+        return name + '%04d' % random.randint(0, 1000)
+    STATE = 1
+else:
+    from pymol.cmd import get_unused_name
+    STATE = -1
+
+
+def _vec_sum(vec_list):
+    # this is the same as
+    # return numpy.array(vec_list).sum(0).tolist()
+    vec = cpv.get_null()
+    for x in vec_list:
+        vec = cpv.add(vec, x)
+    return vec
+
+
+def _mean_and_std(x):
+    # this is the same as
+    # return (numpy.mean(x), numpy.std(x, ddof=1))
+    N = len(x)
+    if N < 2:
+        return (x[0], 0.0)
+    mu = sum(x) / float(N)
+    var = sum([(i - mu) ** 2 for i in x]) / float(N - 1)
+    return (mu, var ** 0.5)
+
+
+def _common_orientation(selection, vec, visualize=1, quiet=0):
+    '''
+    Common part of different helix orientation functions. Does calculate
+    the center of mass and does the visual feedback.
+    '''
+    stored.x = []
+    cmd.iterate_state(STATE, '(%s) and name CA' % (selection),
+                      'stored.x.append([x,y,z])')
+    if len(stored.x) < 2:
+        print('warning: count(CA) < 2')
+        raise CmdException
+    center = cpv.scale(_vec_sum(stored.x), 1. / len(stored.x))
+    if visualize:
+        scale = cpv.distance(stored.x[0], stored.x[-1])
+        visualize_orientation(vec, center, scale, True)
+        cmd.zoom(selection, buffer=2)
+    if not quiet:
+        print('Center: (%.2f, %.2f, %.2f) Direction: (%.2f, %.2f, %.2f)' % tuple(center + vec))
+    return center, vec
+
+
+def visualize_orientation(direction, center=[0, 0, 0], scale=1.0, symmetric=False, color='green', color2='red'):
+    '''
+    Draw an arrow. Helper function for "helix_orientation" etc.
+    '''
+    from pymol import cgo
+    color_list = cmd.get_color_tuple(color)
+    color2_list = cmd.get_color_tuple(color2)
+    if symmetric:
+        scale *= 0.5
+    end = cpv.add(center, cpv.scale(direction, scale))
+    radius = 0.3
+    obj = [cgo.SAUSAGE]
+    obj.extend(center)
+    obj.extend(end)
+    obj.extend([
+        radius,
+        0.8, 0.8, 0.8,
+    ])
+    obj.extend(color_list)
+    if symmetric:
+        start = cpv.sub(center, cpv.scale(direction, scale))
+        obj.append(cgo.SAUSAGE)
+        obj.extend(center)
+        obj.extend(start)
+        obj.extend([
+            radius,
+            0.8, 0.8, 0.8,
+        ])
+        obj.extend(color2_list)
+    coneend = cpv.add(end, cpv.scale(direction, 4.0 * radius))
+    if cmd.get_version()[1] >= 1.2:
+        obj.append(cgo.CONE)
+        obj.extend(end)
+        obj.extend(coneend)
+        obj.extend([
+            radius * 1.75,
+            0.0,
+        ])
+        obj.extend(color_list * 2)
+        obj.extend([
+            1.0, 1.0,  # Caps
+        ])
+    cmd.load_cgo(obj, get_unused_name('oriVec'), zoom=0)
+
+
+def cafit_orientation(selection, visualize=1, quiet=0):
+    '''
+DESCRIPTION
+
+    Get the center and direction of a peptide by least squares
+    linear fit on CA atoms.
+
+USAGE
+
+    cafit_orientation selection [, visualize]
+
+NOTES
+
+    Requires python module "numpy".
+
+SEE ALSO
+
+    helix_orientation
+    '''
+    visualize, quiet = int(visualize), int(quiet)
+    import numpy
+    stored.x = list()
+    cmd.iterate_state(STATE, '(%s) and name CA' % (selection),
+                      'stored.x.append([x,y,z])')
+    x = numpy.array(stored.x)
+    U, s, Vh = numpy.linalg.svd(x - x.mean(0))
+    vec = cpv.normalize(Vh[0])
+    if cpv.dot_product(vec, x[-1] - x[0]) < 0:
+        vec = cpv.negate(vec)
+    return _common_orientation(selection, vec, visualize, quiet)
+
+
+def loop_orientation(selection, visualize=1, quiet=0):
+    '''
+DESCRIPTION
+
+    Get the center and approximate direction of a peptide. Works for any
+    secondary structure.
+    Averages direction of N(i)->C(i) pseudo bonds.
+
+USAGE
+
+    loop_orientation selection [, visualize]
+
+SEE ALSO
+
+    helix_orientation
+    '''
+    visualize, quiet = int(visualize), int(quiet)
+    stored.x = dict()
+    cmd.iterate_state(STATE, '(%s) and name N+C' % (selection),
+                      'stored.x.setdefault(chain + resi, dict())[name] = x,y,z')
+    vec = cpv.get_null()
+    count = 0
+    for x in stored.x.values():
+        if 'C' in x and 'N' in x:
+            vec = cpv.add(vec, cpv.sub(x['C'], x['N']))
+            count += 1
+    if count == 0:
+        print('warning: count == 0')
+        raise CmdException
+    vec = cpv.normalize(vec)
+    return _common_orientation(selection, vec, visualize, quiet)
+
+
+def helix_orientation(selection, visualize=1, sigma_cutoff=1.5, quiet=0):
+    '''
+DESCRIPTION
+
+    Get the center and direction of a helix as vectors. Will only work
+    for helices and gives slightly different results than loop_orientation.
+    Averages direction of C(i)->O(i) bonds.
+
+USAGE
+
+    helix_orientation selection [, visualize [, sigma_cutoff]]
+
+ARGUMENTS
+
+    selection = string: atom selection of helix
+
+    visualize = 0 or 1: show fitted vector as arrow {default: 1}
+
+    sigma_cutoff = float: drop outliers outside
+    (standard_deviation * sigma_cutoff) {default: 1.5}
+
+SEE ALSO
+
+    angle_between_helices, helix_orientation_hbond, loop_orientation, cafit_orientation
+    '''
+    visualize, quiet, sigma_cutoff = int(visualize), int(quiet), float(sigma_cutoff)
+    stored.x = dict()
+    cmd.iterate_state(STATE, '(%s) and name C+O' % (selection),
+                      'stored.x.setdefault(chain + resi, dict())[name] = x,y,z')
+    vec_list = []
+    count = 0
+    for x in stored.x.values():
+        if 'C' in x and 'O' in x:
+            vec_list.append(cpv.sub(x['O'], x['C']))
+            count += 1
+    if count == 0:
+        print('warning: count == 0')
+        raise CmdException
+    vec = _vec_sum(vec_list)
+    if count > 2 and sigma_cutoff > 0:
+        angle_list = [cpv.get_angle(vec, x) for x in vec_list]
+        angle_mu, angle_sigma = _mean_and_std(angle_list)
+        vec_list = [vec_list[i] for i in range(len(vec_list))
+                    if abs(angle_list[i] - angle_mu) < angle_sigma * sigma_cutoff]
+        if not quiet:
+            print('Dropping %d outlier(s)' % (len(angle_list) - len(vec_list)))
+        vec = _vec_sum(vec_list)
+    vec = cpv.normalize(vec)
+    return _common_orientation(selection, vec, visualize, quiet)
+
+
+def helix_orientation_hbond(selection, visualize=1, cutoff=3.5, quiet=0):
+    '''
+DESCRIPTION
+
+    Get the center and direction of a helix as vectors. Will only work
+    for alpha helices and gives slightly different results than
+    helix_orientation. Averages direction of O(i)->N(i+4) hydrogen bonds.
+
+USAGE
+
+    helix_orientation selection [, visualize [, cutoff]]
+
+ARGUMENTS
+
+    cutoff = float: maximal hydrogen bond distance {default: 3.5}
+
+SEE ALSO
+
+    helix_orientation
+    '''
+    visualize, quiet, cutoff = int(visualize), int(quiet), float(cutoff)
+    stored.x = dict()
+    cmd.iterate_state(STATE, '(%s) and name N+O' % (selection),
+                      'stored.x.setdefault(resv, dict())[name] = x,y,z')
+    vec_list = []
+    for resi in stored.x:
+        resi_other = resi + 4
+        if 'O' in stored.x[resi] and resi_other in stored.x:
+            if 'N' in stored.x[resi_other]:
+                vec = cpv.sub(stored.x[resi_other]['N'], stored.x[resi]['O'])
+                if cpv.length(vec) < cutoff:
+                    vec_list.append(vec)
+    if len(vec_list) == 0:
+        print('warning: count == 0')
+        raise CmdException
+    vec = _vec_sum(vec_list)
+    vec = cpv.normalize(vec)
+    return _common_orientation(selection, vec, visualize, quiet)
+
+
+def angle_between_helices(selection1, selection2, method='helix_orientation', visualize=1, quiet=0):
+    '''
+DESCRIPTION
+
+    Calculates the angle between two helices
+
+USAGE
+
+    angle_between_helices selection1, selection2 [, method [, visualize]]
+
+ARGUMENTS
+
+    selection1 = string: atom selection of first helix
+
+    selection2 = string: atom selection of second helix
+
+    method = string: function to calculate orientation {default: helix_orientation}
+             or int: 0: helix_orientation, 1: helix_orientation_hbond,
+                     2: loop_orientation, 3: cafit_orientation
+
+    visualize = 0 or 1: show fitted vector as arrow {default: 1}
+
+SEE ALSO
+
+    helix_orientation, helix_orientation_hbond, loop_orientation, cafit_orientation
+    '''
+    visualize, quiet = int(visualize), int(quiet)
+    methods = {
+        '0': helix_orientation,
+        '1': helix_orientation_hbond,
+        '2': loop_orientation,
+        '3': cafit_orientation,
+    }
+    methods.update([(x.__name__, x) for x in list(methods.values())])
+    try:
+        orientation = methods[str(method)]
+    except KeyError:
+        print('no such method: ' + str(method))
+        raise CmdException
+    if not quiet:
+        print('Using method: ' + orientation.__name__)
+    cen1, dir1 = orientation(selection1, visualize, quiet=1)
+    cen2, dir2 = orientation(selection2, visualize, quiet=1)
+    angle = cpv.get_angle(dir1, dir2)
+    angle_deg = math.degrees(angle)
+    if not quiet:
+        print('Angle: %.2f deg' % (angle_deg))
+    if visualize:
+        cmd.zoom('(%s) or (%s)' % (selection1, selection2), buffer=2)
+    return angle_deg
+
+cmd.extend('helix_orientation', helix_orientation)
+cmd.extend('helix_orientation_hbond', helix_orientation_hbond)
+cmd.extend('loop_orientation', loop_orientation)
+cmd.extend('cafit_orientation', cafit_orientation)
+cmd.extend('angle_between_helices', angle_between_helices)

.pymol/startup/axes.py

@@ -37,8 +37,8 @@ class PutCenterCallback(object):
             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)
+        m = [z, 0, 0, 0, 0, z, 0, 0, 0, 0, z, 0, t[0] / z, t[1] / z, t[2] / z, 1]
+        cmd.set_object_ttt(self.name, m)
 
 def axes(name='axes'):
     '''

.pymol/startup/tmalign.py

@@ -1,257 +0,0 @@
-'''
-(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'],
-})