Difference between revisions of "User:Speleo3"

From PyMOLWiki
Jump to navigation Jump to search
(polarpairs script)
(FreeBSD)
(18 intermediate revisions by the same user not shown)
Line 1: Line 1:
My name is Thomas Holder and I am a bioinformatician at the MPI for Developmental Biology in Tübingen, Germany.
+
My name is Thomas Holder and I work for PyMOL at [http://schrodinger.com Schrödinger].
  
=== Contact ===
+
I was awarded the [http://pymol.org/fellowship Warren L. DeLano Memorial PyMOL Open-Source Fellowship] for 2011-2012.
 +
 
 +
== Contact ==
  
 
* speleo3/users.sourceforge.net
 
* speleo3/users.sourceforge.net
* thomas.holder/tuebingen.mpg.de
+
* thomas.holder/schrodinger.com
  
=== Scripts written by me ===
+
== Scripts written by me ==
  
 
* [[AAindex]]
 
* [[AAindex]]
 
* [[AngleBetweenHelices]]
 
* [[AngleBetweenHelices]]
 +
* [[Extra fit]]
 +
* [[PluginDirectory]]
 +
* [[Pml2py]]
 +
* [[Polarpairs]]
 +
* [[Save settings]]
 +
* [[Show bumps]]
 
* [[Sidechaincenters]]
 
* [[Sidechaincenters]]
 
* [[Spectrumany]]
 
* [[Spectrumany]]
 +
* [[Spectrum states]]
 
* [[Supercell]]
 
* [[Supercell]]
  
=== Scripts Pastebin ===
+
== Scripts Pastebin ==
  
 
Some random scripts with no dedicated PyMOLWiki page.
 
Some random scripts with no dedicated PyMOLWiki page.
 +
 +
Launch interactive python terminal with PyMOL process:
 +
(see also [[Launching From a Script]])
  
 
<source lang="python">
 
<source lang="python">
def extra_fit(selection='(all)', reference=None, method=cmd.align):
+
#!/usr/bin/ipython2.7 -i
     '''
+
 
DESCRIPTION
+
import threading
 +
import pymol._cmd
 +
 
 +
pymol.invocation.parse_args(['pymol', '-qc'])
 +
 
 +
with threading.RLock():
 +
    _COb = pymol._cmd._new(pymol, pymol.invocation.options)
 +
     pymol._cmd._start(_COb, pymol.cmd)
 +
    pymol.cmd._COb = _COb
  
    Like "intra_fit", but for multiple objects instead of
+
from pymol import cmd
    multiple states.
 
    '''
 
    models = cmd.get_object_list(selection)
 
    if reference is None:
 
        reference = models[0]
 
        models = models[1:]
 
    elif reference in models:
 
        models.remove(reference)
 
    if cmd.is_string(method):
 
        method = eval(method)
 
    for model in models:
 
        print model, method(model, reference)
 
cmd.extend('extra_fit', extra_fit)
 
 
</source>
 
</source>
  
 
<source lang="python">
 
<source lang="python">
def mse2met(selection='all', quiet=1):
+
#!/usr/bin/python2.6 -i
    '''
+
 
DESCRIPTION
+
import sys, os
 +
 
 +
# autocompletion
 +
import readline
 +
import rlcompleter
 +
readline.parse_and_bind('tab: complete')
 +
 
 +
# pymol environment
 +
moddir='/opt/pymol-svn/modules'
 +
sys.path.insert(0, moddir)
 +
os.putenv('PYMOL_PATH', os.path.join(moddir, 'pymol/pymol_path'))
 +
 
 +
# pymol launching
 +
import pymol
 +
pymol.pymol_argv = ['pymol','-qc'] + sys.argv[1:]
 +
pymol.finish_launching()
 +
cmd = pymol.cmd
 +
</source>
 +
 
 +
Build FREEMOL (see also [[MovieSchool 6]])
 +
 
 +
<source lang="bash">
 +
#!/bin/bash -e
 +
 
 +
src=/tmp
 +
prefix=/opt/pymol-git
 +
export FREEMOL=$prefix/freemol
 +
 
 +
freemoltrunk=$src/freemol-trunk
 +
if [[ ! -e $freemoltrunk ]]; then
 +
    svn co svn://bioinformatics.org/svnroot/freemol/trunk $freemoltrunk
 +
fi
 +
 
 +
cd $freemoltrunk
 +
 
 +
sed -i 's/vdwtype\[11\]/vdwtype[14]/' src/mengine/src/field.h
 +
 
 +
for name in mpeg_encode mengine apbs pdb2pqr; do
 +
    (cd src/$name && ./configure && make && make install)
 +
done
 +
 
 +
cp -na freemol/libpy/freemol $prefix/modules/
 +
 
 +
ln -sfT $FREEMOL $prefix/modules/pymol/pymol_path/freemol
 +
</source>
 +
 
 +
Download all PyMOL scripts from Robert L. Campbell's website:
 +
 
 +
<source lang="bash">
 +
wget -r -np -nd --level=1 -A .py \
 +
    http://pldserver1.biochem.queensu.ca/~rlc/work/pymol/
 +
</source>
 +
 
 +
Render movie from PNG files (save as <code>png2mpeg1.sh</code>):
 +
 
 +
<source lang="bash">
 +
#!/bin/bash
 +
 
 +
set -e
 +
 
 +
usage="usage: $(basename $0) [-w width] [-f fps] [-b vbitrate] <indir> <outfile.mpeg>"
 +
width=""
 +
fps=25
 +
vbitrate=16000
 +
 
 +
args="$(getopt w:f:b:h "$@")" || args="-h"
 +
set -- $args
 +
 
 +
while [[ $# > 0 ]]; do
 +
    case "$1" in
 +
        --) shift; break ;;
 +
        -w) width=$2; shift 2 ;;
 +
        -f) fps=$2; shift 2 ;;
 +
        -b) vbitrate=$2; shift 2 ;;
 +
        -h) echo $usage; exit 1 ;;
 +
        *) echo "argument error: $1"; exit 1 ;;
 +
    esac
 +
done
 +
 
 +
if [[ $# > 2 ]]; then
 +
    echo "too many arguments: $3 ..."
 +
    echo $usage
 +
    exit 1
 +
fi
  
    Mutate selenomethionine to methionine
+
indir="$1"
    '''
+
outfile="$2"
    quiet = int(quiet)
 
    x = cmd.alter('(%s) and MSE/SE' % selection, 'name="SD";elem="S"')
 
    cmd.alter('(%s) and MSE/' % selection, 'resn="MET";type="ATOM"')
 
    if not quiet:
 
        print 'Altered %d MSE residues to MET' % (x)
 
    cmd.sort()
 
cmd.extend('mse2met', mse2met)
 
  
def remove_alt(selection='all', keep='A', quiet=1):
+
if [[ -z "$indir" ]]; then
     '''
+
    echo "error: indir missing"
DESCRIPTION
+
    echo $usage
 +
     exit 2
 +
fi
  
     Remove alternative location atoms.
+
if [[ -z "$outfile" ]]; then
 +
     echo "error: outfile missing"
 +
    echo $usage
 +
    exit 3
 +
fi
  
ARGUMENTS
+
MENCODER="mencoder -quiet"
 +
MPEG1ARGS="-mf type=png:fps=$fps -ovc lavc -forceidx -noskip \
 +
    -of rawvideo -mpegopts format=mpeg1 \
 +
    -lavcopts vcodec=mpeg1video:vbitrate=$vbitrate:vhq:trell:keyint=25"
  
     selection = string: atom selection
+
if [[ -n "$width" ]]; then
 +
     MPEG1ARGS="-zoom -xy $width -sws 9 $MPEG1ARGS"
 +
fi
  
    keep = string: AltLoc to keep {default: A}
+
pattern="mf://$indir/*.png"
    '''
+
$MENCODER "$pattern" $MPEG1ARGS:vpass=1 -o /dev/null
    cmd.remove('(%s) and not alt +%s' % (selection, keep), quiet=int(quiet))
+
$MENCODER "$pattern" $MPEG1ARGS:vpass=2 -o "$outfile"
    cmd.alter(selection, 'alt=""')
 
    cmd.sort()
 
cmd.extend('remove_alt', remove_alt)
 
 
</source>
 
</source>
 +
 +
=== Export movie with transparent background ===
  
 
<source lang="python">
 
<source lang="python">
def cbm(selection='(all)'):
+
# make transparent pngs
    '''
+
set opaque_background, off
DESCRIPTION
+
set ray_trace_frames
 +
mpng foo
 +
 
 +
# create movie file (use codec "qtrle" or "png")
 +
system ffmpeg -i foo%04d.png -vcodec qtrle foo.mov
 +
 
 +
# clean up
 +
system rm -f foo????.png
 +
</source>
 +
 
 +
=== load_mtz_cctbx: Load MTZ files with a [[CCTBX]] wrapper (3 files) ===
 +
 
 +
1) ~/bin/mtz2ccp4.sh
  
    Color by molecule
+
<source lang="bash">
    '''
+
#!/bin/bash
    col = 2
+
export PATH=/opt/ccp4/ccp4-6.5/bin:$PATH
    for model in cmd.get_object_list(selection):
+
exec cctbx.python ~/bin/mtz2ccp4.py "$@"
        cmd.color(col, '%s and (%s)' % (model, selection))
 
        col += 1
 
cmd.extend('cbm', cbm)
 
 
</source>
 
</source>
 +
 +
2) ~/bin/mtz2ccp4.py
  
 
<source lang="python">
 
<source lang="python">
def dev2b(selection='name CA'):
+
#!/opt/ccp4/ccp4-6.5/bin/cctbx.python
 +
 
 +
import os
 +
import sys
 +
import tempfile
 +
 
 +
def mtz2ccp4maps(filename, prefix='map'):
 +
    '''
 +
Creates a temporary directory and dumps all maps from the given MTZ file
 +
into this directory as CCP4 maps files. Returns the path of the temporary
 +
directory.
 
     '''
 
     '''
DESCRIPTION
+
    from iotbx.reflection_file_reader import any_reflection_file
 +
 
 +
    hkl_in = any_reflection_file(file_name=filename)
 +
 
 +
    temp_dir = tempfile.mkdtemp()
 +
 
 +
    for i_map, array in enumerate(hkl_in.as_miller_arrays()):
 +
        if array.is_complex_array():
 +
            fft_map = array.fft_map(resolution_factor=0.25).apply_sigma_scaling()
 +
            map_filename = os.path.join(temp_dir,
 +
                    prefix + '_' + '_'.join(array.info().labels) + '.ccp4')
 +
            fft_map.as_ccp4_map(file_name=map_filename)
 +
 
 +
    return temp_dir
  
    Determine the RMSD per residue for a multi-state object and assign b-factor
+
# print the name of the temporary directory to standard output
    '''
+
print mtz2ccp4maps(*sys.argv[1:])
    import numpy
 
    stored.x = dict()
 
    stored.b = dict()
 
    for state in range(1, cmd.count_states()+1):
 
        cmd.iterate_state(state, selection, 'stored.x.setdefault((model,segi,chain,resi,name), []).append([x,y,z])')
 
    for key, coord_list in stored.x.iteritems():
 
        b_sq = numpy.array(coord_list).var(0).mean() # var over states, mean over x,y,z
 
        stored.b[key] = numpy.sqrt(b_sq) * 10.0
 
    cmd.alter(selection, 'b = stored.b[model,segi,chain,resi,name]')
 
cmd.extend('dev2b', dev2b)
 
 
</source>
 
</source>
 +
 +
3) ~/.pymolrc.py
  
 
<source lang="python">
 
<source lang="python">
def select_extendbyss(selection, name=None, quiet=0):
+
@cmd.extend
 +
def load_mtz_cctbx(filename, prefix=''):
 
     '''
 
     '''
 
DESCRIPTION
 
DESCRIPTION
  
     Extend selection by connected secondary structure elements.
+
     Load all maps from an MTZ file, using the mtz2ccp4.sh wrapper which
 +
    uses iotbx (cctbx).
 +
    '''
 +
    import subprocess
 +
    import glob
 +
    import shutil
  
ARGUMENTS
+
    if not prefix:
 +
        prefix = os.path.basename(filename).rpartition('.')[0]
 +
 
 +
    outdir = subprocess.Popen([os.path.expanduser('~/bin/mtz2ccp4.sh'),
 +
        filename, prefix], stdout=subprocess.PIPE).stdout.readlines()[0].strip()
  
     selection = string: selection-expression
+
     for mapfilename in glob.glob(os.path.join(outdir, '*.ccp4')):
 +
        cmd.load(mapfilename)
  
     name = string: create a named atom selection if not None {default: None}
+
     shutil.rmtree(outdir)
    '''
 
    def in_intervals(i, intervals):
 
        for interval in intervals:
 
            if interval[0] <= i and i <= interval[1]:
 
                return True
 
        return False
 
    quiet = int(quiet)
 
    stored.x = set()
 
    # only iterate over CAs since for example PyMOL's dss command just
 
    # assignes ss to CAs.
 
    cmd.iterate('bycalpha (%s)' % (selection),
 
            'stored.x.add((model,segi,chain,resv,ss))')
 
    elements = dict()
 
    for model,segi,chain,resv,ss in stored.x:
 
        key = (model,segi,chain,ss)
 
        elements.setdefault(key, [])
 
        if in_intervals(resv, elements[key]):
 
            continue
 
        stored.b = set()
 
        cmd.iterate('/%s/%s/%s//CA and ss "%s"' % key, 'stored.b.add(resv)')
 
        resv_min = resv
 
        resv_max = resv
 
        while (resv_min - 1) in stored.b:
 
            resv_min -= 1
 
        while (resv_max + 1) in stored.b:
 
            resv_max += 1
 
        elements[key].append((resv_min, resv_max))
 
    sele_list = []
 
    ss_names = {'S': 'Strand', 'H': 'Helix', '': 'Loop'}
 
    for key in elements:
 
        model,segi,chain,ss = key
 
        for resv_min,resv_max in elements[key]:
 
            sele = '/%s/%s/%s/%d-%d' % (model, segi, chain, resv_min, resv_max)
 
            sele_list.append(sele)
 
            if not quiet:
 
                print ss_names.get(ss, ss), sele
 
    sele = ' or '.join(sele_list)
 
    if name is not None:
 
        cmd.select(name, sele)
 
    if not quiet:
 
        print 'Selection:', sele
 
    return sele
 
cmd.extend('select_extendbyss', select_extendbyss)
 
 
</source>
 
</source>
 +
 +
=== ccmutate ===
  
 
<source lang="python">
 
<source lang="python">
def diff(sele1, sele2, byres=1, name=None, operator='in', quiet=0):
+
@cmd.extend
 +
def ccmutate(code, selection='??sele|?pk1', sculpt=1):
 
     '''
 
     '''
 
DESCRIPTION
 
DESCRIPTION
  
     Difference between two molecules
+
     Mutate selected residue.
  
 
ARGUMENTS
 
ARGUMENTS
  
     sele1 = string: atom selection
+
     code = str: 3-letter PDBeChem chemical component identifier
 +
 
 +
    selection = str: single residue selection {default: pk1 or sele}
  
     sele2 = string: atom selection
+
     sculpt = 0/1: try to adopt conformation of replaced sidechain, followed
 +
    by relaxation using sculpting {default: 1}
  
    byres = 0/1: report residues, not atoms (does not affect selection)
+
EXAMPLE
    {default: 1}
 
  
     operator = in/like/align: operator to match atoms {default: in}
+
     fetch 1ubq, async=0
 +
    ccmutate 0HG, resi 24
  
 
SEE ALSO
 
SEE ALSO
  
     symdiff
+
     fetch ..., type=cc
 +
    wizard mutagenesis
 
     '''
 
     '''
     byres, quiet = int(byres), int(quiet)
+
     code = code.upper()
     if name is None:
+
 
        name = cmd.get_unused_name('diff')
+
     tmp_sele = cmd.get_unused_name('_sele')
     if operator == 'align':
+
     tmp_frag = cmd.get_unused_name('_frag')
        alnobj = cmd.get_unused_name('__aln')
+
    tmp_Nnbr = cmd.get_unused_name('_Nnbr')
        cmd.align(sele1, sele2, cycles=0, transform=0, object=alnobj)
+
    tmp_back = cmd.get_unused_name('_back')
        sele = '(%s) and not %s' % (sele1, alnobj)
+
    tmp_tmpl = cmd.get_unused_name('_tmpl')
        cmd.select(name, sele)
+
    tmp_sc_o = cmd.get_unused_name('_sc_o')
        cmd.delete(alnobj)
+
    tmp_sc_n = cmd.get_unused_name('_sc_n')
     else:
+
 
         sele = '(%s) and not ((%s) %s (%s))' % (sele1, sele1, operator, sele2)
+
     try:
         cmd.select(name, sele)
+
         cmd.select(tmp_sele, 'byres (' + selection + ')', 0)
    if not quiet:
+
 
        if byres:
+
         # check input selection
             seleiter = 'byca ' + name
+
        if cmd.count_atoms('name CA & ?' + tmp_sele) != 1:
            expr = 'print "/%s/%s/%s/%s`%s" % (model,segi,chain,resn,resi)'
+
             raise pymol.CmdException('selection must include exactly one residue')
        else:
+
        if cmd.count_atoms('name N+CA+C & ?' + tmp_sele) != 3:
            seleiter = name
+
            raise pymol.CmdException("selected residue doesn't have N+CA+C atoms")
            expr = 'print "/%s/%s/%s/%s`%s/%s" % (model,segi,chain,resn,resi,name)'
 
        cmd.iterate(seleiter, expr)
 
    return name
 
  
def symdiff(sele1, sele2, byres=1, name=None, operator='in', quiet=0):
+
        # PDBeChem fragment
    '''
+
        cmd.fetch(code, tmp_frag, type='cc', zoom=0)
DESCRIPTION
+
 
 +
        # check if fragment is amino acid
 +
        if cmd.count_atoms('name N+CA+C & ?' + tmp_frag) != 3:
 +
            raise pymol.CmdException("residue '%s' doesn't have N+CA+C atoms" % (code))
  
    Symmetric difference between two molecules
+
        # only keep hydrogens if target also has hydrogens
 +
        if cmd.count_atoms('hydro & ?' + tmp_sele) == 0:
 +
            cmd.remove('hydro & ?' + tmp_frag)
  
SEE ALSO
+
        # update residue name for old residue
 +
        cmd.alter(tmp_sele, 'resn = ' + repr(code))
  
    diff
+
        # superpose fragment on backbone
    '''
+
         cmd.align(
    byres, quiet = int(byres), int(quiet)
+
                'name N+CA+C & ?' + tmp_frag,
    if name is None:
+
                'name N+CA+C & ?' + tmp_sele)
         name = cmd.get_unused_name('symdiff')
 
    tmpname = cmd.get_unused_name('__tmp')
 
    diff(sele1, sele2, byres, name, operator, quiet)
 
    diff(sele2, sele1, byres, tmpname, operator, quiet)
 
    cmd.select(name, tmpname, merge=1)
 
    cmd.delete(tmpname)
 
    return name
 
  
cmd.extend('symdiff', symdiff)
+
        # extra N bonds, like in PRO
cmd.extend('diff', diff)
+
        cmd.select(tmp_Nnbr, 'neighbor (name N & ?' + tmp_frag + ')', 0)
</source>
 
  
<source lang="python">
+
        # backbone selection
def polarpairs(sel1, sel2, cutoff=4.0, angle=63.0, name='', quiet=1):
+
        cmd.select(tmp_back, 'name CA+C+O+N+OXT', 0)
    '''
+
        cmd.select(tmp_back, 'hydro & neighbor ?' + tmp_back, 0, merge=1)
ARGUMENTS
 
  
    sel1, sel2 = string: atom selections
+
        # remove complementary atoms
 +
        cmd.remove(          '?' + tmp_frag + ' & ' + tmp_back)
 +
        cmd.extract(tmp_tmpl, '?' + tmp_sele + ' & !' + tmp_back, zoom=0)
  
    cutoff = float: distance cutoff
+
        if cmd.count_atoms(tmp_frag):
 +
            # attach new sidechain
 +
            cmd.fuse('name CB & ?' + tmp_frag, 'name CA & ?' + tmp_sele, mode=1, move=0)
 +
            cmd.unpick()
  
    angle = float: h-bond angle cutoff in degrees. If angle="default", take
+
            # new atom selections
    "h_bond_max_angle" setting. If angle=0, do not detect h-bonding.
+
            cmd.select(tmp_sc_n, '(byres ?' + tmp_sele + ') & !?' + tmp_sele, 0)
 +
            cmd.select(tmp_sc_o, '?' + tmp_sc_n + ' like ?' + tmp_tmpl, 0)
  
SEE ALSO
+
            # extra N bonds, like in PRO
 +
            if cmd.count_atoms(tmp_Nnbr):
 +
                cmd.bond('?' + tmp_sc_n + ' like ?' + tmp_Nnbr, 'name N & ?' + tmp_sele)
  
    cmd.find_pairs, cmd.distance
+
            # adopt old conformation, if possible
    '''
+
            if int(sculpt):
    cutoff = float(cutoff)
+
                model = cmd.get_object_list('?' + tmp_sele)[0]
    quiet = int(quiet)
+
                cmd.protect(model)
    if angle == 'default':
+
                cmd.deprotect('?' + tmp_sc_n + ' & !?' + tmp_sc_o)
        angle = cmd.get('h_bond_max_angle', sel1)
+
                cmd.sculpt_activate(model)
    angle = float(angle)
+
                if cmd.count_atoms(tmp_sc_o):
    mode = 1 if angle > 0 else 0
+
                    cmd.update(tmp_sc_o, tmp_tmpl)
    x = cmd.find_pairs('(%s) and donors' % sel1, '(%s) and acceptors' % sel2,
+
                    cmd.set('sculpt_field_mask', 63) # local geom + vdw
            cutoff=cutoff, mode=mode, angle=angle) + \
+
                    cmd.sculpt_iterate(model, cycles=100)
        cmd.find_pairs('(%s) and acceptors' % sel1, '(%s) and donors' % sel2,
+
                    cmd.deprotect(tmp_sc_o)
            cutoff=cutoff, mode=mode, angle=angle)
+
                cmd.set('sculpt_field_mask', 0xff) # all
    x = sorted(set(x))
+
                cmd.sculpt_iterate(model, cycles=200)
    if not quiet:
+
                cmd.set('sculpt_field_mask', 31) # local geom
        print 'Settings: cutoff=%.1fangstrom angle=%.1fdegree' % (cutoff, angle)
+
                cmd.sculpt_iterate(model, cycles=200)
        print 'Found %d polar contacts' % (len(x))
 
    if len(name) > 0:
 
        for p in x:
 
            cmd.distance(name, '(%s`%s)' % p[0], '(%s`%s)' % p[1])
 
    return x
 
  
cmd.extend('polarpairs', polarpairs)
+
    finally:
 +
        cmd.delete(tmp_sele)
 +
        cmd.delete(tmp_frag)
 +
        cmd.delete(tmp_Nnbr)
 +
        cmd.delete(tmp_back)
 +
        cmd.delete(tmp_tmpl)
 +
        cmd.delete(tmp_sc_o)
 +
        cmd.delete(tmp_sc_n)
 
</source>
 
</source>
  
Launch interactive python terminal with PyMOL process:
+
=== CTRL-L ligand zoom ===
  
<source lang="bash">
+
<source lang="python">
#!/bin/bash
+
@cmd.set_key('CTRL-L')
moddir=/opt/pymol-svn/lib/python2.6/site-packages
+
def ligand_zoom():
export PYTHONPATH="$moddir:$PYTHONPATH"
+
    global _current_ligand
export PYMOL_PATH="$moddir/pymol/pymol_path"
+
    s = {'ligand_set': set()}
export PYMOL_DATA="$PYMOL_PATH/data"
+
    if cmd.iterate('organic', 'ligand_set.add((model,segi,chain,resi))',
exec python2.6 -i -c "
+
            space=s) < 1:
import readline
+
        return
import rlcompleter
+
    ligands = sorted(s["ligand_set"])
readline.parse_and_bind('tab: complete')
+
    try:
import pymol
+
        i = ligands.index(_current_ligand)
pymol.pymol_argv = ['pymol','-qc']
+
    except (ValueError, NameError):
pymol.finish_launching()
+
        i = -1
cmd = pymol.cmd
+
    i = (i + 1) % len(ligands)
"
+
    _current_ligand = ligands[i]
 +
    # use "do" for feedback
 +
    cmd.do('zoom /%s/%s/%s & resi %s, animate=1, buffer=2' % ligands[i])
 
</source>
 
</source>
 +
 +
=== Compile on FreeBSD ===
 +
 +
pkg upgrade
 +
pkg install subversion py27-Pmw glew freeglut png freetype2 libxml2 msgpack
 +
python2 setup.py install --prefix=$HOME/opt/pymol-svn

Revision as of 14:20, 16 December 2017

My name is Thomas Holder and I work for PyMOL at Schrödinger.

I was awarded the Warren L. DeLano Memorial PyMOL Open-Source Fellowship for 2011-2012.

Contact

  • speleo3/users.sourceforge.net
  • thomas.holder/schrodinger.com

Scripts written by me

Scripts Pastebin

Some random scripts with no dedicated PyMOLWiki page.

Launch interactive python terminal with PyMOL process: (see also Launching From a Script)

#!/usr/bin/ipython2.7 -i

import threading
import pymol._cmd

pymol.invocation.parse_args(['pymol', '-qc'])

with threading.RLock():
    _COb = pymol._cmd._new(pymol, pymol.invocation.options)
    pymol._cmd._start(_COb, pymol.cmd)
    pymol.cmd._COb = _COb

from pymol import cmd
#!/usr/bin/python2.6 -i

import sys, os

# autocompletion
import readline
import rlcompleter
readline.parse_and_bind('tab: complete')

# pymol environment
moddir='/opt/pymol-svn/modules'
sys.path.insert(0, moddir)
os.putenv('PYMOL_PATH', os.path.join(moddir, 'pymol/pymol_path'))

# pymol launching
import pymol
pymol.pymol_argv = ['pymol','-qc'] + sys.argv[1:]
pymol.finish_launching()
cmd = pymol.cmd

Build FREEMOL (see also MovieSchool 6)

#!/bin/bash -e

src=/tmp
prefix=/opt/pymol-git
export FREEMOL=$prefix/freemol

freemoltrunk=$src/freemol-trunk
if [[ ! -e $freemoltrunk ]]; then
    svn co svn://bioinformatics.org/svnroot/freemol/trunk $freemoltrunk
fi

cd $freemoltrunk

sed -i 's/vdwtype\[11\]/vdwtype[14]/' src/mengine/src/field.h

for name in mpeg_encode mengine apbs pdb2pqr; do
    (cd src/$name && ./configure && make && make install)
done

cp -na freemol/libpy/freemol $prefix/modules/

ln -sfT $FREEMOL $prefix/modules/pymol/pymol_path/freemol

Download all PyMOL scripts from Robert L. Campbell's website:

wget -r -np -nd --level=1 -A .py \
    http://pldserver1.biochem.queensu.ca/~rlc/work/pymol/

Render movie from PNG files (save as png2mpeg1.sh):

#!/bin/bash

set -e

usage="usage: $(basename $0) [-w width] [-f fps] [-b vbitrate] <indir> <outfile.mpeg>"
width=""
fps=25
vbitrate=16000

args="$(getopt w:f:b:h "$@")" || args="-h"
set -- $args

while [[ $# > 0 ]]; do
    case "$1" in
        --) shift; break ;;
        -w) width=$2; shift 2 ;;
        -f) fps=$2; shift 2 ;;
        -b) vbitrate=$2; shift 2 ;;
        -h) echo $usage; exit 1 ;;
        *) echo "argument error: $1"; exit 1 ;;
    esac
done

if [[ $# > 2 ]]; then
    echo "too many arguments: $3 ..."
    echo $usage
    exit 1
fi

indir="$1"
outfile="$2"

if [[ -z "$indir" ]]; then
    echo "error: indir missing"
    echo $usage
    exit 2
fi

if [[ -z "$outfile" ]]; then
    echo "error: outfile missing"
    echo $usage
    exit 3
fi

MENCODER="mencoder -quiet"
MPEG1ARGS="-mf type=png:fps=$fps -ovc lavc -forceidx -noskip \
    -of rawvideo -mpegopts format=mpeg1 \
    -lavcopts vcodec=mpeg1video:vbitrate=$vbitrate:vhq:trell:keyint=25"

if [[ -n "$width" ]]; then
    MPEG1ARGS="-zoom -xy $width -sws 9 $MPEG1ARGS"
fi

pattern="mf://$indir/*.png"
$MENCODER "$pattern" $MPEG1ARGS:vpass=1 -o /dev/null
$MENCODER "$pattern" $MPEG1ARGS:vpass=2 -o "$outfile"

Export movie with transparent background

# make transparent pngs
set opaque_background, off
set ray_trace_frames
mpng foo

# create movie file (use codec "qtrle" or "png")
system ffmpeg -i foo%04d.png -vcodec qtrle foo.mov 

# clean up
system rm -f foo????.png

load_mtz_cctbx: Load MTZ files with a CCTBX wrapper (3 files)

1) ~/bin/mtz2ccp4.sh

#!/bin/bash
export PATH=/opt/ccp4/ccp4-6.5/bin:$PATH
exec cctbx.python ~/bin/mtz2ccp4.py "$@"

2) ~/bin/mtz2ccp4.py

#!/opt/ccp4/ccp4-6.5/bin/cctbx.python

import os
import sys
import tempfile

def mtz2ccp4maps(filename, prefix='map'):
    '''
Creates a temporary directory and dumps all maps from the given MTZ file
into this directory as CCP4 maps files. Returns the path of the temporary
directory.
    '''
    from iotbx.reflection_file_reader import any_reflection_file

    hkl_in = any_reflection_file(file_name=filename)

    temp_dir = tempfile.mkdtemp()

    for i_map, array in enumerate(hkl_in.as_miller_arrays()):
        if array.is_complex_array():
            fft_map = array.fft_map(resolution_factor=0.25).apply_sigma_scaling()
            map_filename = os.path.join(temp_dir,
                    prefix + '_' + '_'.join(array.info().labels) + '.ccp4')
            fft_map.as_ccp4_map(file_name=map_filename)

    return temp_dir

# print the name of the temporary directory to standard output
print mtz2ccp4maps(*sys.argv[1:])

3) ~/.pymolrc.py

@cmd.extend
def load_mtz_cctbx(filename, prefix=''):
    '''
DESCRIPTION

    Load all maps from an MTZ file, using the mtz2ccp4.sh wrapper which
    uses iotbx (cctbx).
    '''
    import subprocess
    import glob
    import shutil

    if not prefix:
        prefix = os.path.basename(filename).rpartition('.')[0]

    outdir = subprocess.Popen([os.path.expanduser('~/bin/mtz2ccp4.sh'),
        filename, prefix], stdout=subprocess.PIPE).stdout.readlines()[0].strip()

    for mapfilename in glob.glob(os.path.join(outdir, '*.ccp4')):
        cmd.load(mapfilename)

    shutil.rmtree(outdir)

ccmutate

@cmd.extend
def ccmutate(code, selection='??sele|?pk1', sculpt=1):
    '''
DESCRIPTION

    Mutate selected residue.

ARGUMENTS

    code = str: 3-letter PDBeChem chemical component identifier

    selection = str: single residue selection {default: pk1 or sele}

    sculpt = 0/1: try to adopt conformation of replaced sidechain, followed
    by relaxation using sculpting {default: 1}

EXAMPLE

    fetch 1ubq, async=0
    ccmutate 0HG, resi 24

SEE ALSO

    fetch ..., type=cc
    wizard mutagenesis
    '''
    code = code.upper()

    tmp_sele = cmd.get_unused_name('_sele')
    tmp_frag = cmd.get_unused_name('_frag')
    tmp_Nnbr = cmd.get_unused_name('_Nnbr')
    tmp_back = cmd.get_unused_name('_back')
    tmp_tmpl = cmd.get_unused_name('_tmpl')
    tmp_sc_o = cmd.get_unused_name('_sc_o')
    tmp_sc_n = cmd.get_unused_name('_sc_n')

    try:
        cmd.select(tmp_sele, 'byres (' + selection + ')', 0)

        # check input selection
        if cmd.count_atoms('name CA & ?' + tmp_sele) != 1:
            raise pymol.CmdException('selection must include exactly one residue')
        if cmd.count_atoms('name N+CA+C & ?' + tmp_sele) != 3:
            raise pymol.CmdException("selected residue doesn't have N+CA+C atoms")

        # PDBeChem fragment
        cmd.fetch(code, tmp_frag, type='cc', zoom=0)

        # check if fragment is amino acid
        if cmd.count_atoms('name N+CA+C & ?' + tmp_frag) != 3:
            raise pymol.CmdException("residue '%s' doesn't have N+CA+C atoms" % (code))

        # only keep hydrogens if target also has hydrogens
        if cmd.count_atoms('hydro & ?' + tmp_sele) == 0:
            cmd.remove('hydro & ?' + tmp_frag)

        # update residue name for old residue
        cmd.alter(tmp_sele, 'resn = ' + repr(code))

        # superpose fragment on backbone
        cmd.align(
                'name N+CA+C & ?' + tmp_frag,
                'name N+CA+C & ?' + tmp_sele)

        # extra N bonds, like in PRO
        cmd.select(tmp_Nnbr, 'neighbor (name N & ?' + tmp_frag + ')', 0)

        # backbone selection
        cmd.select(tmp_back, 'name CA+C+O+N+OXT', 0)
        cmd.select(tmp_back, 'hydro & neighbor ?' + tmp_back, 0, merge=1)

        # remove complementary atoms
        cmd.remove(           '?' + tmp_frag + ' & ' + tmp_back)
        cmd.extract(tmp_tmpl, '?' + tmp_sele + ' & !' + tmp_back, zoom=0)

        if cmd.count_atoms(tmp_frag):
            # attach new sidechain
            cmd.fuse('name CB & ?' + tmp_frag, 'name CA & ?' + tmp_sele, mode=1, move=0)
            cmd.unpick()

            # new atom selections
            cmd.select(tmp_sc_n, '(byres ?' + tmp_sele + ') & !?' + tmp_sele, 0)
            cmd.select(tmp_sc_o, '?' + tmp_sc_n + ' like ?' + tmp_tmpl, 0)

            # extra N bonds, like in PRO
            if cmd.count_atoms(tmp_Nnbr):
                cmd.bond('?' + tmp_sc_n + ' like ?' + tmp_Nnbr, 'name N & ?' + tmp_sele)

            # adopt old conformation, if possible
            if int(sculpt):
                model = cmd.get_object_list('?' + tmp_sele)[0]
                cmd.protect(model)
                cmd.deprotect('?' + tmp_sc_n + ' & !?' + tmp_sc_o)
                cmd.sculpt_activate(model)
                if cmd.count_atoms(tmp_sc_o):
                    cmd.update(tmp_sc_o, tmp_tmpl)
                    cmd.set('sculpt_field_mask', 63) # local geom + vdw
                    cmd.sculpt_iterate(model, cycles=100)
                    cmd.deprotect(tmp_sc_o)
                cmd.set('sculpt_field_mask', 0xff) # all
                cmd.sculpt_iterate(model, cycles=200)
                cmd.set('sculpt_field_mask', 31) # local geom
                cmd.sculpt_iterate(model, cycles=200)

    finally:
        cmd.delete(tmp_sele)
        cmd.delete(tmp_frag)
        cmd.delete(tmp_Nnbr)
        cmd.delete(tmp_back)
        cmd.delete(tmp_tmpl)
        cmd.delete(tmp_sc_o)
        cmd.delete(tmp_sc_n)

CTRL-L ligand zoom

@cmd.set_key('CTRL-L')
def ligand_zoom():
    global _current_ligand
    s = {'ligand_set': set()}
    if cmd.iterate('organic', 'ligand_set.add((model,segi,chain,resi))',
            space=s) < 1:
        return
    ligands = sorted(s["ligand_set"])
    try:
        i = ligands.index(_current_ligand)
    except (ValueError, NameError):
        i = -1
    i = (i + 1) % len(ligands)
    _current_ligand = ligands[i]
    # use "do" for feedback
    cmd.do('zoom /%s/%s/%s & resi %s, animate=1, buffer=2' % ligands[i])

Compile on FreeBSD

pkg upgrade
pkg install subversion py27-Pmw glew freeglut png freetype2 libxml2 msgpack
python2 setup.py install --prefix=$HOME/opt/pymol-svn