# -*- coding: utf-8 -*-
"""This module defines miscellaneous functions dealing with protein data."""
import numpy as np
from prody.atomic import Atomic, Atom, AtomGroup, Selection, HierView
from prody.utilities import openFile, showFigure
from prody import SETTINGS
__all__ = ['showProtein', 'writePQR', ]
[docs]def writePQR(filename, atoms):
"""Write *atoms* in PQR format to a file with name *filename*. Only
current coordinate set is written. Returns *filename* upon success. If
*filename* ends with :file:`.gz`, a compressed file will be written."""
if not isinstance(atoms, Atomic):
raise TypeError('atoms does not have a valid type')
if isinstance(atoms, Atom):
atoms = Selection(atoms.getAtomGroup(), [atoms.getIndex()],
atoms.getACSIndex(),
'index ' + str(atoms.getIndex()))
stream = openFile(filename, 'w')
n_atoms = atoms.numAtoms()
atomnames = atoms.getNames()
if atomnames is None:
raise RuntimeError('atom names are not set')
for i, an in enumerate(atomnames):
lenan = len(an)
if lenan < 4:
atomnames[i] = ' ' + an
elif lenan > 4:
atomnames[i] = an[:4]
s_or_u = np.array(['a']).dtype.char
resnames = atoms._getResnames()
if resnames is None:
resnames = ['UNK'] * n_atoms
resnums = atoms._getResnums()
if resnums is None:
resnums = np.ones(n_atoms, int)
chainids = atoms._getChids()
if chainids is None:
chainids = np.zeros(n_atoms, s_or_u + '1')
charges = atoms._getCharges()
if charges is None:
charges = np.zeros(n_atoms, float)
radii = atoms._getRadii()
if radii is None:
radii = np.zeros(n_atoms, float)
icodes = atoms._getIcodes()
if icodes is None:
icodes = np.zeros(n_atoms, s_or_u + '1')
hetero = ['ATOM'] * n_atoms
heteroflags = atoms._getFlags('hetatm')
if heteroflags is None:
heteroflags = atoms._getFlags('hetero')
if heteroflags is not None:
hetero = np.array(hetero, s_or_u + '6')
hetero[heteroflags] = 'HETATM'
altlocs = atoms._getAltlocs()
if altlocs is None:
altlocs = np.zeros(n_atoms, s_or_u + '1')
format = ('{0:6s}{1:5d} {2:4s}{3:1s}' +
'{4:4s}{5:1s}{6:4d}{7:1s} ' +
'{8:8.3f}{9:8.3f}{10:8.3f}' +
'{11:8.4f}{12:7.4f}\n').format
coords = atoms._getCoords()
write = stream.write
for i, xyz in enumerate(coords):
write(format(hetero[i], i+1, atomnames[i], altlocs[i],
resnames[i], chainids[i], int(resnums[i]),
icodes[i], xyz[0], xyz[1], xyz[2], charges[i], radii[i]))
write('TER\nEND')
stream.close()
return filename
[docs]def showProtein(*atoms, **kwargs):
"""Show protein representation using :meth:`~mpl_toolkits.mplot3d.Axes3D`.
This function is designed for generating a quick view of the contents of a
:class:`~.AtomGroup` or :class:`~.Selection`.
Protein atoms matching ``"calpha"`` selection are displayed using solid
lines by picking a random and unique color per chain. Line with can
be adjusted using *lw* argument, e.g. ``lw=12``. Default width is 4.
Chain colors can be overwritten using chain identifier as in ``A='green'``.
Water molecule oxygen atoms are represented by red colored circles. Color
can be changed using *water* keyword argument, e.g. ``water='aqua'``.
Water marker and size can be changed using *wmarker* and *wsize* keywords,
defaults values are ``wmarker='.', wsize=6``.
Hetero atoms matching ``"hetero and noh"`` selection are represented by
circles and unique colors are picked at random on a per residue basis.
Colors can be customized using residue name as in ``NAH='purple'``. Note
that this will color all distinct residues with the same name in the same
color. Hetero atom marker and size can be changed using *hmarker* and
*hsize* keywords, default values are ``hmarker='o', hsize=6``.
ProDy will set the size of axis so the representation is not distorted when
the shape of figure window is close to a square. Colors are picked at
random, except for water oxygens which will always be colored red."""
alist = atoms
for atoms in alist:
if not isinstance(atoms, Atomic):
raise TypeError('atoms must be an Atomic instance')
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
cf = plt.gcf()
show = None
for child in cf.get_children():
if isinstance(child, Axes3D):
show = child
break
if show is None:
show = Axes3D(cf)
from matplotlib import colors
cnames = dict(colors.cnames)
wcolor = kwargs.get('water', 'red').lower()
avoid = np.array(colors.hex2color(cnames.pop(wcolor, cnames.pop('red'))))
for cn, val in cnames.items(): # PY3K: OK
clr = np.array(colors.hex2color(val))
if clr.sum() > 2.4:
cnames.pop(cn)
elif np.abs(avoid - clr).sum() <= 0.6:
cnames.pop(cn)
cnames = list(cnames)
import random
random.shuffle(cnames)
min_ = list()
max_ = list()
for atoms in alist:
if isinstance(atoms, AtomGroup):
title = atoms.getTitle()
else:
title = atoms.getAtomGroup().getTitle()
calpha = atoms.select('calpha')
if calpha:
for ch in HierView(calpha, chain=True):
xyz = ch._getCoords()
chid = ch.getChid()
show.plot(xyz[:, 0], xyz[:, 1], xyz[:, 2],
label=title + '_' + chid,
color=kwargs.get(chid, cnames.pop()).lower(),
lw=kwargs.get('lw', 4))
water = atoms.select('water and noh')
if water:
xyz = atoms.select('water')._getCoords()
show.plot(xyz[:, 0], xyz[:, 1], xyz[:, 2], label=title + '_water',
color=wcolor,
ls='None', marker=kwargs.get('wmarker', '.'),
ms=kwargs.get('wsize', 6))
hetero = atoms.select('not protein and not nucleic and not water')
if hetero:
for res in HierView(hetero).iterResidues():
xyz = res._getCoords()
resname = res.getResname()
resnum = str(res.getResnum())
chid = res.getChid()
show.plot(xyz[:, 0], xyz[:, 1], xyz[:, 2], ls='None',
color=kwargs.get(resname, cnames.pop()).lower(),
label=title + '_' + chid + '_' + resname + resnum,
marker=kwargs.get('hmarker', 'o'),
ms=kwargs.get('hsize', 6))
xyz = atoms._getCoords()
min_.append(xyz.min(0))
max_.append(xyz.max(0))
show.set_xlabel('x')
show.set_ylabel('y')
show.set_zlabel('z')
min_ = np.array(min_).min(0)
max_ = np.array(max_).max(0)
center = (max_ + min_) / 2
half = (max_ - min_).max() / 2
show.set_xlim3d(center[0]-half, center[0]+half)
show.set_ylim3d(center[1]-half, center[1]+half)
show.set_zlim3d(center[2]-half, center[2]+half)
if kwargs.get('legend', False):
show.legend(prop={'size': 10})
if SETTINGS['auto_show']:
showFigure()
return show