~ubuntu-branches/ubuntu/natty/python-cogent/natty

« back to all changes in this revision

Viewing changes to doc/cookbook/structural_data.rst

  • Committer: Bazaar Package Importer
  • Author(s): Steffen Moeller
  • Date: 2010-12-04 22:30:35 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20101204223035-j11kinhcrrdgg2p2
Tags: 1.5-1
* Bumped standard to 3.9.1, no changes required.
* New upstream version.
  - major additions to Cookbook
  - added AlleleFreqs attribute to ensembl Variation objects.
  - added getGeneByStableId method to genome objects.
  - added Introns attribute to Transcript objects and an Intron class.
  - added Mann-Whitney test and a Monte-Carlo version
  - exploratory and confirmatory period estimation techniques (suitable for
    symbolic and continuous data)
  - Information theoretic measures (AIC and BIC) added
  - drawing of trees with collapsed nodes
  - progress display indicator support for terminal and GUI apps
  - added parser for illumina HiSeq2000 and GAiix sequence files as 
    cogent.parse.illumina_sequence.MinimalIlluminaSequenceParser.
  - added parser to FASTQ files, one of the output options for illumina's
    workflow, also added cookbook demo.
  - added functionality for parsing of SFF files without the Roche tools in
    cogent.parse.binary_sff
  - thousand fold performance improvement to nmds
  - >10-fold performance improvements to some Table operations

Show diffs side-by-side

added added

removed removed

Lines of Context:
1
 
***************
2
1
Structural data
3
 
***************
 
2
---------------
4
3
 
5
4
.. sectionauthor:: Kristian Rother, Patrick Yannul
6
5
 
7
6
Protein structures
8
 
==================
 
7
^^^^^^^^^^^^^^^^^^
 
8
 
 
9
Reading Protein structures
 
10
""""""""""""""""""""""""""
 
11
 
 
12
Retrieve a structure from PDB
 
13
+++++++++++++++++++++++++++++
 
14
 
 
15
.. doctest::
 
16
 
 
17
    >>> from cogent.db.pdb import Pdb
 
18
    >>> p = Pdb()
 
19
    >>> pdb_file = p['4tsv']
 
20
    >>> pdb = pdb_file.read()
 
21
    >>> len(pdb)
 
22
    135027
 
23
 
 
24
This example will retrieve the structure as a PDB file string.
 
25
 
 
26
Parse a PDB file
 
27
++++++++++++++++
 
28
 
 
29
.. doctest::
 
30
 
 
31
    >>> from cogent.parse.pdb import PDBParser
 
32
    >>> struc = PDBParser(open('data/4TSV.pdb'))
 
33
    >>> struc
 
34
    <Structure id=4TSV>
 
35
 
 
36
Parse a PDB entry directly from the web
 
37
+++++++++++++++++++++++++++++++++++++++
 
38
 
 
39
.. doctest::
 
40
 
 
41
    >>> from cogent.parse.pdb import PDBParser
 
42
    >>> struc = PDBParser(p['4tsv'])
 
43
 
 
44
Accessing PDB header information
 
45
++++++++++++++++++++++++++++++++
 
46
 
 
47
.. doctest::
 
48
 
 
49
    >>> struc.header['id']
 
50
    '4TSV'
 
51
    >>> struc.header['resolution']
 
52
    '1.80'
 
53
    >>> struc.header['r_free']
 
54
    '0.262'
 
55
    >>> struc.header['space_group']
 
56
    'H 3'
 
57
 
 
58
Navigating structure objects
 
59
""""""""""""""""""""""""""""
 
60
 
 
61
What does a structure object contain?
 
62
+++++++++++++++++++++++++++++++++++++
 
63
 
 
64
A ``cogent.parse.pdb.Structure`` object as returned by ``PDBParser`` contains a tree-like hierarchy of ``Entity`` objects. They are organized such that ``Structures`` that contain ``Models`` that contain ``Chains`` that contain ``Residues`` that in turn contain ``Atoms``. You can read more about the entity model on [URL of Marcins example page].
 
65
 
 
66
How to access a model from a structure
 
67
++++++++++++++++++++++++++++++++++++++
 
68
 
 
69
To get the first model out of a structure:
 
70
 
 
71
.. doctest::
 
72
 
 
73
    >>> model = struc[(0,)]
 
74
    >>> model
 
75
    <Model id=0>
 
76
 
 
77
The key contains the model number as a tuple.
 
78
 
 
79
How to access a chain from a model?
 
80
+++++++++++++++++++++++++++++++++++
 
81
 
 
82
To get a particular chain:
 
83
 
 
84
.. doctest::
 
85
 
 
86
    >>> chain = model[('A',)]
 
87
    >>> chain
 
88
    <Chain id=A>
 
89
 
 
90
How to access a residue from a chain?
 
91
+++++++++++++++++++++++++++++++++++++
 
92
 
 
93
To get a particular residue:
 
94
 
 
95
.. doctest::
 
96
 
 
97
    >>> resi = chain[('ILE', 154, ' '),]
 
98
    >>> resi
 
99
    <Residue ILE resseq=154 icode= >
 
100
 
 
101
What properties does a residue have?
 
102
++++++++++++++++++++++++++++++++++++
 
103
 
 
104
.. doctest::
 
105
 
 
106
    >>> resi.res_id
 
107
    154
 
108
    >>> resi.name
 
109
    'ILE'
 
110
    >>> resi.h_flag
 
111
    ' '
 
112
    >>> resi.seg_id
 
113
    '    '
 
114
 
 
115
 
 
116
Access an atom from a residue
 
117
+++++++++++++++++++++++++++++
 
118
 
 
119
To get a particular atom:
 
120
 
 
121
.. doctest::
 
122
 
 
123
    >>> atom = resi[("N", ' '),]
 
124
    >>> atom
 
125
    <Atom ('N', ' ')>
 
126
 
 
127
Properties of an atom
 
128
+++++++++++++++++++++
 
129
 
 
130
.. doctest::
 
131
 
 
132
    >>> atom.name
 
133
    ' N  '
 
134
    >>> atom.element
 
135
    ' N'
 
136
    >>> atom.coords
 
137
    array([ 142.986,   36.523,    6.838])
 
138
    >>> atom.bfactor
 
139
    13.35
 
140
    >>> atom.occupancy
 
141
    1.0
 
142
 
 
143
If a model/chain/residue/atom does not exist
 
144
++++++++++++++++++++++++++++++++++++++++++++
 
145
 
 
146
You will get a ``KeyError``.
 
147
 
 
148
Is there something special about heteroatoms to consider?
 
149
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
150
 
 
151
Yes, they have the ``h_flag`` attribute set in residues.
 
152
 
 
153
How are Altlocs/insertion codes represented?
 
154
++++++++++++++++++++++++++++++++++++++++++++
 
155
 
 
156
Both are part of the residue/atom ID.
 
157
 
 
158
Useful methods to access Structure objects
 
159
""""""""""""""""""""""""""""""""""""""""""
 
160
 
 
161
How to access all atoms, residues etc via a dictionary
 
162
++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
163
 
 
164
The ``table`` property of a structure returns a two-dimensional dictionary containing all atoms. The keys are 1) the entity level (any of 'A','R','C','M') and 2) the combined IDs of ``Structure``, ``Model``, ``Chain``, ``Residue``, ``Atom`` as a tuple.
 
165
 
 
166
.. doctest::
 
167
 
 
168
    >>> struc.table['A'][('4TSV', 0, 'A', ('HIS', 73, ' '), ('O', ' '))]
 
169
    <Atom ('O', ' ')>
 
170
 
 
171
Calculate the center of mass of a model or chain
 
172
++++++++++++++++++++++++++++++++++++++++++++++++
 
173
 
 
174
.. NEEDS TO BE CHECKED WITH MARCIN
 
175
 
 
176
.. doctest::
 
177
 
 
178
    >>> model.coords
 
179
    array([ 146.66615752,   35.08673503,   -3.60735847])
 
180
    >>> chain.coords
 
181
    array([ 146.66615752,   35.08673503,   -3.60735847])
 
182
 
 
183
 
 
184
How to get a list of all residues in a chain?
 
185
+++++++++++++++++++++++++++++++++++++++++++++
 
186
 
 
187
.. doctest::
 
188
 
 
189
    >>> chain.values()[0]
 
190
    <Residue ILE resseq=154 icode= >
 
191
 
 
192
How to get a list of all atoms in a chain?
 
193
++++++++++++++++++++++++++++++++++++++++++
 
194
 
 
195
.. doctest::
 
196
 
 
197
    >>> resi.values()[0]
 
198
    <Atom ('N', ' ')>
 
199
 
 
200
Constructing structures
 
201
"""""""""""""""""""""""
 
202
 
 
203
How to create a new entity?
 
204
+++++++++++++++++++++++++++
 
205
 
 
206
``Structure``/``Model``/``Chain``/``Residue``/``Atom`` objects can be created as follows:
 
207
 
 
208
.. doctest::
 
209
 
 
210
    >>> from cogent.core.entity import Structure,Model,Chain,Residue,Atom
 
211
    >>> from numpy import array
 
212
    >>> s = Structure('my_struc')
 
213
    >>> m = Model((0),)
 
214
    >>> c = Chain(('A'),)
 
215
    >>> r = Residue(('ALA', 1, ' ',),False,' ')
 
216
    >>> a = Atom(('C  ',' ',), 'C', 1, array([0.0,0.0,0.0]), 1.0, 0.0, 'C')
 
217
 
 
218
How to add entities to each other?
 
219
++++++++++++++++++++++++++++++++++
 
220
 
 
221
.. doctest::
 
222
 
 
223
    >>> s.addChild(m)
 
224
    >>> m.addChild(c)
 
225
    >>> c.addChild(r)
 
226
    >>> r.addChild(a)
 
227
    >>> s.setTable(force=True)
 
228
    >>> s.table
 
229
    {'A': {('my_struc', 0, 'A', ('ALA', 1, ' '), ('C  ', ' ')): <Atom ('C  ', ' ')>}, 'C': {('my_struc', 0, 'A'): <Chain id=A>}, 'R': {('my_struc', 0, 'A', ('ALA', 1, ' ')): <Residue ALA resseq=1 icode= >}, 'M': {('my_struc', 0): <Model id=0>}}
 
230
 
 
231
How to remove a residue from a chain?
 
232
+++++++++++++++++++++++++++++++++++++
 
233
 
 
234
.. doctest::
 
235
 
 
236
    >>> c.delChild(r.id)
 
237
    >>> s.table
 
238
    {'A': {('my_struc', 0, 'A', ('ALA', 1, ' '), ...
 
239
 
 
240
Geometrical analyses
 
241
""""""""""""""""""""
9
242
 
10
243
Calculating euclidean distances between atoms
11
 
---------------------------------------------
12
 
 
13
 
.. doctest::
14
 
    
 
244
+++++++++++++++++++++++++++++++++++++++++++++
 
245
 
 
246
.. doctest::
 
247
 
 
248
    >>> from cogent.maths.geometry import distance
 
249
    >>> atom1 = resi[('N', ' '),]
 
250
    >>> atom2 = resi[('CA', ' '),]
 
251
    >>> distance(atom1.coords, atom2.coords)
 
252
    1.4691967192993618
 
253
 
 
254
Calculating euclidean distances between coordinates
 
255
+++++++++++++++++++++++++++++++++++++++++++++++++++
 
256
 
 
257
.. doctest::
 
258
 
15
259
    >>> from numpy import array
16
260
    >>> from cogent.maths.geometry import distance
17
261
    >>> a1 = array([1.0, 2.0, 3.0])
18
262
    >>> a2 = array([1.0, 4.0, 9.0])
19
263
    >>> distance(a1,a2)
20
264
    6.324...
 
265
 
 
266
Calculating flat angles from atoms
 
267
++++++++++++++++++++++++++++++++++
 
268
 
 
269
.. doctest::
 
270
 
 
271
    >>> from cogent.struct.dihedral import angle
 
272
    >>> atom3 = resi[('C', ' '),]
 
273
    >>> a12 = atom2.coords-atom1.coords
 
274
    >>> a23 = atom2.coords-atom3.coords
 
275
    >>> angle(a12,a23)
 
276
    1.856818...
 
277
 
 
278
Calculates the angle in radians.
 
279
 
 
280
Calculating flat angles from coordinates
 
281
++++++++++++++++++++++++++++++++++++++++
 
282
 
 
283
.. doctest::
 
284
 
 
285
    >>> from cogent.struct.dihedral import angle
 
286
    >>> a1 = array([0.0, 0.0, 1.0])
 
287
    >>> a2 = array([0.0, 0.0, 0.0])
 
288
    >>> a3 = array([0.0, 1.0, 0.0])
 
289
    >>> a12 = a2-a1
 
290
    >>> a23 = a2-a3
 
291
    >>> angle(a12,a23)
 
292
    1.5707963267948966
 
293
 
 
294
Calculates the angle in radians.
 
295
 
 
296
Calculating dihedral angles from atoms
 
297
++++++++++++++++++++++++++++++++++++++
 
298
 
 
299
.. doctest::
 
300
 
 
301
    >>> from cogent.struct.dihedral import dihedral
 
302
    >>> atom4 = resi[('CG1', ' '),]
 
303
    >>> dihedral(atom1.coords,atom2.coords,atom3.coords, atom4.coords)
 
304
    259.49277688244217
 
305
 
 
306
Calculates the torsion in degrees.
 
307
 
 
308
Calculating dihedral angles from coordinates
 
309
++++++++++++++++++++++++++++++++++++++++++++
 
310
 
 
311
.. doctest::
 
312
 
 
313
    >>> from cogent.struct.dihedral import dihedral
 
314
    >>> a1 = array([0.0, 0.0, 1.0])
 
315
    >>> a2 = array([0.0, 0.0, 0.0])
 
316
    >>> a3 = array([0.0, 1.0, 0.0])
 
317
    >>> a4 = array([1.0, 1.0, 0.0])
 
318
    >>> dihedral(a1,a2,a3,a4)
 
319
    90.0
 
320
 
 
321
Calculates the torsion in degrees.
 
322
 
 
323
Other stuff
 
324
"""""""""""
 
325
 
 
326
How to count the atoms in a structure?
 
327
++++++++++++++++++++++++++++++++++++++
 
328
 
 
329
.. doctest::
 
330
 
 
331
    >>> len(struc.table['A'].values())
 
332
    1187
 
333
 
 
334
How to iterate over chains in canonical PDB order?
 
335
++++++++++++++++++++++++++++++++++++++++++++++++++
 
336
 
 
337
In PDB, the chain with space as ID comes last, the
 
338
others in alphabetical order.
 
339
 
 
340
.. doctest::
 
341
 
 
342
    >>> for chain in model.sortedvalues():
 
343
    ...     print chain
 
344
    <Chain id=A>
 
345
 
 
346
How to iterate over chains in alphabetical order?
 
347
+++++++++++++++++++++++++++++++++++++++++++++++++
 
348
 
 
349
If you want the chains in purely alphabetical order:
 
350
 
 
351
.. KR 2 ROB: Is this what you requested or is the above example enough?
 
352
 
 
353
.. doctest::
 
354
 
 
355
    >>> keys = model.keys()
 
356
    >>> keys.sort()
 
357
    >>> for chain in [model[id] for id in keys]:
 
358
    ...     print chain
 
359
    <Chain id=A>
 
360
 
 
361
How to iterate over all residues in a chain?
 
362
++++++++++++++++++++++++++++++++++++++++++++
 
363
 
 
364
.. doctest::
 
365
 
 
366
    >>> residues = [resi for resi in chain.values()]
 
367
    >>> len(residues)
 
368
    218
 
369
 
 
370
How to remove all water molecules from a structure
 
371
++++++++++++++++++++++++++++++++++++++++++++++++++
 
372
 
 
373
.. doctest::
 
374
 
 
375
    >>> water = [r for r in struc.table['R'].values() if r.name == 'H_HOH']
 
376
    >>> for resi in water:
 
377
    ...     resi.parent.delChild(resi.id)
 
378
    >>> struc.setTable(force=True)
 
379
    >>> len(struc.table['A'].values())
 
380
    1117
 
381
    >>> residues = [resi for resi in chain.values()]
 
382
    >>> len(residues)
 
383
    148
 
384