-
Notifications
You must be signed in to change notification settings - Fork 0
/
rest.py
101 lines (83 loc) · 2.93 KB
/
rest.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
# Author: Gael Varoquaux <gael.varoquaux@normalesup.org>
# Copyright (c) 2008-2015, Enthought, Inc.
# License: BSD Style.
# The pdb code for the protein.
protein_code = '2q09'
# Retrieve the file from the protein database #################################
import os
if not os.path.exists('pdb%s.ent.gz' % protein_code):
# Download the data
try:
from urllib import urlopen
except ImportError:
from urllib.request import urlopen
print('Downloading protein data, please wait')
opener = urlopen(
'ftp://ftp.wwpdb.org/pub/pdb/data/structures/divided/pdb/q0/pdb%s.ent.gz'
% protein_code)
open('pdb%s.ent.gz' % protein_code, 'wb').write(opener.read())
# Parse the pdb file ##########################################################
import gzip
infile = gzip.GzipFile('pdb%s.ent.gz' % protein_code, 'rb')
# A graph represented by a dictionary associating nodes with keys
# (numbers), and edges (pairs of node keys).
nodes = dict()
edges = list()
atoms = set()
# Build the graph from the PDB information
last_atom_label = None
last_chain_label = None
for line in infile:
line = line.split()
if line[0] in ('ATOM', 'HETATM'):
nodes[line[1]] = (line[2], line[6], line[7], line[8])
atoms.add(line[2])
chain_label = line[5]
if chain_label == last_chain_label:
edges.append((line[1], last_atom_label))
last_atom_label = line[1]
last_chain_label = chain_label
elif line[0] == 'CONECT':
for start, stop in zip(line[1:-1], line[2:]):
edges.append((start, stop))
atoms = list(atoms)
atoms.sort()
atoms = dict(zip(atoms, range(len(atoms))))
# Turn the graph into 3D positions, and a connection list.
labels = dict()
x = list()
y = list()
z = list()
scalars = list()
for index, label in enumerate(nodes):
labels[label] = index
this_scalar, this_x, this_y, this_z = nodes[label]
scalars.append(atoms[this_scalar])
x.append(float(this_x))
y.append(float(this_y))
z.append(float(this_z))
connections = list()
for start, stop in edges:
connections.append((labels[start], labels[stop]))
import numpy as np
x = np.array(x)
y = np.array(y)
z = np.array(z)
scalars = np.array(scalars)
# Visualize the data ##########################################################
from mayavi import mlab
mlab.figure(1, bgcolor=(0, 0, 0))
mlab.clf()
pts = mlab.points3d(x, y, z, 1.5 * scalars.max() - scalars,
scale_factor=0.015, resolution=10)
pts.mlab_source.dataset.lines = np.array(connections)
# Use a tube fiter to plot tubes on the link, varying the radius with the
# scalar value
tube = mlab.pipeline.tube(pts, tube_radius=0.15)
tube.filter.radius_factor = 1.
tube.filter.vary_radius = 'vary_radius_by_scalar'
mlab.pipeline.surface(tube, color=(0.8, 0.8, 0))
# Visualize the local atomic density
mlab.pipeline.volume(mlab.pipeline.gaussian_splatter(pts))
mlab.view(49, 31.5, 52.8, (4.2, 37.3, 20.6))
mlab.show()