1 """
2 3D and secondary structure APIs.
3
4 This module defines some of the most fundamental abstractions in the library:
5 L{Structure}, L{Chain}, L{Residue} and L{Atom}. Instances of these objects may
6 exist independently and that is perfectly fine, but usually they are part of a
7 Composite aggregation. The root node in this Composite is a L{Structure} (or
8 L{Ensemble}). L{Structure}s are composed of L{Chain}s, and each L{Chain} is a
9 collection of L{Residue}s. The leaf node is L{Atom}.
10
11 All of these objects implement the base L{AbstractEntity} interface. Therefore,
12 every node in the Composite can be transformed:
13
14 >>> r, t = [rotation matrix], [translation vector]
15 >>> entity.transform(r, t)
16
17 and it knows its immediate children:
18
19 >>> entity.items
20 <iterator> # over all immediate child entities
21
22 If you want to traverse the complete Composite tree, starting at arbitrary level,
23 and down to the lowest level, use one of the L{CompositeEntityIterator}s. Or just
24 call L{AbstractEntity.components}:
25
26 >>> entity.components()
27 <iterator> # over all descendants, of any type, at any level
28 >>> entity.components(klass=Residue)
29 <iterator> # over all Residue descendants
30
31 Some of the inner objects in this hierarchy behave just like dictionaries
32 (but are not):
33
34 >>> structure.chains['A'] # access chain A by ID
35 <Chain A: Protein>
36 >>> structure['A'] # the same
37 <Chain A: Protein>
38 >>> residue.atoms['CS']
39 <Atom: CA> # access an atom by its name
40 >>> residue.atoms['CS']
41 <Atom: CA> # the same
42
43 Others behave like list collections:
44
45 >>> chain.residues[10] # 1-based access to the residues in the chain
46 <ProteinResidue [10]: PRO 10>
47 >>> chain[10] # 0-based, list-like access
48 <ProteinResidue [11]: GLY 11>
49
50 Step-wise building of L{Ensemble}s, L{Chain}s and L{Residue}s is supported through
51 a number of C{append} methods, for example:
52
53 >>> residue = ProteinResidue(401, ProteinAlphabet.ALA)
54 >>> s.chains['A'].residues.append(residue)
55
56 See L{EnsembleModelsCollection}, L{StructureChainsTable}, L{ChainResiduesCollection}
57 and L{ResidueAtomsTable} for more details.
58
59 Some other objects in this module of potential interest are the self-explanatory
60 L{SecondaryStructure} and L{TorsionAngles}.
61 """
62
63 import os
64 import re
65 import copy
66 import math
67 import numpy
68
69 import csb.io
70 import csb.core
71 import csb.numeric
72 import csb.bio.utils
73
74 from abc import ABCMeta, abstractmethod, abstractproperty
75
76 from csb.bio.sequence import SequenceTypes, SequenceAlphabets, AlignmentTypes
80 """
81 Torsion angle unit types
82 """
83 Degrees='deg'; Radians='rad'
84
91
105
109
112
115
118
121
124
127
130
133
136
139
142
145
148
150 """
151 Base class for all protein structure entities.
152
153 This class defines uniform interface of all entities (e.g. L{Structure},
154 L{Chain}, L{Residue}) according to the Composite pattern.
155 """
156
157 __metaclass__ = ABCMeta
158
159 @abstractproperty
161 """
162 Iterator over all immediate children of the entity
163 @rtype: iterator of L{AbstractEntity}
164 """
165 pass
166
168 """
169 Return an iterator over all descendants of the entity.
170
171 @param klass: return entities of the specified L{AbstractEntity} subclass
172 only. If None, traverse the hierarchy down to the lowest level.
173 @param klass: class
174 """
175 for entity in CompositeEntityIterator.create(self, klass):
176 if klass is None or isinstance(entity, klass):
177 yield entity
178
188
190 """
191 Extract the coordinates of the specified kind(s) of atoms and return
192 them as a list.
193
194 @param what: a list of atom kinds, e.g. ['N', 'CA', 'C']
195 @type what: list or None
196
197 @return: a list of lists, each internal list corresponding to the coordinates
198 of a 3D vector
199 @rtype: list
200
201 @raise Broken3DStructureError: if a specific atom kind cannot be retrieved from a residue
202 """
203 coords = [ ]
204
205 for residue in self.components(klass=Residue):
206 for atom_kind in (what or residue.atoms):
207 try:
208 coords.append(residue.atoms[atom_kind].vector)
209 except csb.core.ItemNotFoundError:
210 if skip:
211 continue
212 raise Broken3DStructureError('Could not retrieve {0} atom from the structure'.format(atom_kind))
213
214 return numpy.array(coords)
215
217 """
218 Iterates over composite L{AbstractEntity} hierarchies.
219
220 @param node: root entity to traverse
221 @type node: L{AbstractEntity}
222 """
223
225
226 if not isinstance(node, AbstractEntity):
227 raise TypeError(node)
228
229 self._node = node
230 self._stack = csb.core.Stack()
231
232 self._inspect(node)
233
236
239
241
242 while True:
243 if self._stack.empty():
244 raise StopIteration()
245
246 try:
247 current = self._stack.peek()
248 node = next(current)
249 self._inspect(node)
250 return node
251
252 except StopIteration:
253 self._stack.pop()
254
256 """
257 Push C{node}'s children to the stack.
258 """
259 self._stack.push(node.items)
260
261 @staticmethod
263 """
264 Create a new composite iterator.
265
266 @param leaf: if not None, return a L{ConfinedEntityIterator}
267 @type leaf: class
268 @rtype: L{CompositeEntityIterator}
269 """
270 if leaf is None:
271 return CompositeEntityIterator(node)
272 else:
273 return ConfinedEntityIterator(node, leaf)
274
276 """
277 Iterates over composite L{AbstractEntity} hierarchies, but terminates
278 the traversal of a branch once a specific node type is encountered.
279
280 @param node: root entity to traverse
281 @type node: L{AbstractEntity}
282 @param leaf: traverse the hierarchy down to the specified L{AbstractEntity}
283 @type leaf: class
284 """
292
294
295 if not isinstance(node, self._leaf):
296 self._stack.push(node.items)
297
298 -class Ensemble(csb.core.AbstractNIContainer, AbstractEntity):
299 """
300 Represents an ensemble of multiple L{Structure} models.
301 Provides a list-like access to these models:
302
303 >>> ensemble[0]
304 <Structure Model 1: accn, x chains>
305 >>> ensemble.models[1]
306 <Structure Model 1: accn, x chains>
307 """
308
311
314
315 @property
318
319 @property
321 """
322 Access Ensembles's models by model ID
323 @rtype: L{EnsembleModelsCollection}
324 """
325 return self._models
326
327 @property
329 return iter(self._models)
330
331 @property
333 """
334 The first L{Structure} in the ensemble (if available)
335 @rtype: L{Structure} or None
336 """
337 if len(self._models) > 0:
338 return self[0]
339 return None
340
341 - def to_pdb(self, output_file=None):
371
373
378
393
394 @property
397
398
399 -class Structure(csb.core.AbstractNIContainer, AbstractEntity):
400 """
401 Represents a single model of a PDB 3-Dimensional molecular structure.
402 Provides access to the L{Chain} objects, contained in the model:
403
404 >>> structure['A']
405 <Chain A: Protein>
406 >>> structure.chains['A']
407 <Chain A: Protein>
408 >>> structure.items
409 <iterator of Chain-s>
410
411 @param accession: accession number of the structure
412 @type accession: str
413 """
422
424 return "<Structure Model {0.model_id}: {0.accession}, {1} chains>".format(self, self.chains.length)
425
426 @property
429
430 @property
432 """
433 Access chains by their chain identifiers
434 @rtype: L{StructureChainsTable}
435 """
436 return self._chains
437
438 @property
440 for chain in self._chains:
441 yield self._chains[chain]
442
443 @property
445 """
446 The first L{Chain} in the structure (if available)
447 @rtype: L{Chain} or None
448 """
449 if len(self._chains) > 0:
450 return next(self.items)
451 return None
452
453 @property
455 """
456 Accession number
457 @rtype: str
458 """
459 return self._accession
460 @accession.setter
467
468 @property
470 """
471 Model ID
472 @rtype: int
473 """
474 return self._model_id
475 @model_id.setter
477 self._model_id = value
478
479 @property
481 """
482 Resolution in Angstroms
483 """
484 return self._resolution
485 @resolution.setter
490
492 """
493 Dump the structure in FASTA format.
494
495 @return: FASTA-formatted string with all chains in the structure
496 @rtype: str
497
498 @deprecated: this method will be removed soon. Use
499 L{csb.bio.sequence.ChainSequence.create} instead
500 """
501 fasta = []
502
503 for chain in self.items:
504
505 if chain.length > 0:
506 fasta.append('>{0}'.format(chain.header))
507 fasta.append(chain.sequence)
508
509 return os.linesep.join(fasta)
510
511 - def to_pdb(self, output_file=None):
535
537
538 - def __init__(self, structure=None, chains=None):
545
547 if len(self) > 0:
548 return "<StructureChains: {0}>".format(', '.join(self))
549 else:
550 return "<StructureChains: empty>"
551
552 @property
555
557 """
558 Add a new Chain to the structure.
559
560 @param chain: the new chain to be appended
561 @type chain: L{Chain}
562
563 @raise DuplicateChainIDError: if a chain with same ID is already defined
564 @raise InvalidOperation: if the chain is already associated with a structure
565 """
566
567 if chain._structure and chain._structure is not self.__container:
568 raise InvalidOperation('This chain is already part of another structure')
569 if chain.id in self:
570 raise DuplicateChainIDError('A chain with ID {0} is already defined'.format(chain.id))
571
572 super(StructureChainsTable, self).append(chain.id, chain)
573
574 if self.__container:
575 chain._accession = self.__container.accession
576 chain._structure = self.__container
577
579 """
580 Remove a chain from the structure.
581
582 @param id: ID of the chain to be detached
583 @type id: str
584 """
585 chain = self[id]
586 self._remove(id)
587 chain._structure = None
588
600
601 -class Chain(csb.core.AbstractNIContainer, AbstractEntity):
602 """
603 Represents a polymeric chain. Provides list-like and rank-based access to
604 the residues in the chain:
605
606 >>> chain[0]
607 <ProteinResidue [1]: SER None>
608 >>> chain.residues[1]
609 <ProteinResidue [1]: SER None>
610
611 You can also access residues by their PDB sequence number:
612
613 >>> chain.find(sequence_number=5, insertion_code='A')
614 <ProteinResidue [1]: SER 5A>
615
616 @param chain_id: ID of the new chain
617 @type chain_id: str
618 @param type: sequence type (a member of the L{SequenceTypes} enum)
619 @type type: L{csb.core.EnumItem}
620 @param name: name of the chain
621 @type name: str
622 @param residues: initialization list of L{Residue}-s
623 @type residues: list
624 @param accession: accession number of the chain
625 @type accession: str
626 @param molecule_id: MOL ID of the chain, if part of a polymer
627
628 """
631
632 self._id = str(chain_id).strip()
633 self._accession = None
634 self._type = None
635 self._residues = ChainResiduesCollection(self, residues)
636 self._secondary_structure = None
637 self._molecule_id = molecule_id
638 self._torsion_computed = False
639 self._name = str(name).strip()
640
641 self._structure = None
642
643 self.type = type
644 if accession is not None:
645 self.accession = accession
646
647 @staticmethod
665
666 @property
668 return self._residues
669
671 return "<Chain {0.id}: {0.type!r}>".format(self)
672
674 return self._residues.length
675
676 @property
678 """
679 Chain's ID
680 @rtype: str
681 """
682 return self._id
683 @id.setter
685 if not isinstance(id, csb.core.string):
686 raise ValueError(id)
687 id = id.strip()
688 if self._structure:
689 self._structure.chains._update_chain_id(self, id)
690 self._id = id
691
692 @property
694 """
695 Accession number
696 @rtype: str
697 """
698 return self._accession
699 @accession.setter
706
707 @property
709 """
710 Chain type - any member of L{SequenceTypes}
711 @rtype: enum item
712 """
713 return self._type
714 @type.setter
715 - def type(self, type):
719
720 @property
722 """
723 Rank-based access to Chain's L{Residue}s
724 @rtype: L{ChainResiduesCollection}
725 """
726 return self._residues
727
728 @property
730 return iter(self._residues)
731
732 @property
751
752 @property
754 """
755 True if C{Chain.compute_torsion} had been invoked
756 @rtype: bool
757 """
758 return self._torsion_computed
759
760 @property
762 """
763 Number of residues
764 @rtype: int
765 """
766 return self._residues.length
767
768 @property
769 - def entry_id(self):
770 """
771 Accession number + chain ID
772 @rtype: str
773 """
774 if self._accession and self._id:
775 return self._accession + self._id
776 else:
777 return None
778
779 @property
781 """
782 Chain name
783 @rtype: str
784 """
785 return self._name
786 @name.setter
787 - def name(self, value):
791
792 @property
794 """
795 PDB MOL ID of this chain
796 @rtype: int
797 """
798 return self._molecule_id
799 @molecule_id.setter
801 self._molecule_id = value
802
803 @property
805 """
806 FASTA header in PDB format
807 @rtype: str
808 """
809 header = "{0._accession}_{0._id} mol:{1} length:{0.length} {0.name}"
810 return header.format(self, str(self.type).lower())
811
812 @property
828
829 @property
831 """
832 Sequence alphabet corresponding to the current chain type
833 @rtype: L{csb.core.enum}
834 """
835 return SequenceAlphabets.get(self.type)
836
837 @property
839 """
840 Secondary structure (if available)
841 @rtype: L{SecondaryStructure}
842 """
843 return self._secondary_structure
844 @secondary_structure.setter
852
863
864 - def subregion(self, start, end, clone=False):
865 """
866 Extract a subchain defined by [start, end]. If clone is True, this
867 is a deep copy of the chain. Otherwise same as:
868
869 >>> chain.residues[start : end + 1]
870
871 but coordinates are checked and a Chain instance is returned.
872
873 @param start: start position of the sub-region
874 @type start: int
875 @param end: end position
876 @type end: int
877 @param clone: if True, a deep copy of the sub-region is returned,
878 otherwise - a shallow one
879 @type clone: bool
880
881
882 @return: a new chain, made from the residues of the extracted region
883 @rtype: L{Chain}
884
885 @raise IndexError: if start/end positions are out of range
886 """
887 if start < self.residues.start_index or start > self.residues.last_index:
888 raise IndexError('The start position is out of range {0.start_index} .. {0.last_index}'.format(self.residues))
889 if end < self.residues.start_index or end > self.residues.last_index:
890 raise IndexError('The end position is out of range {0.start_index} .. {0.last_index}'.format(self.residues))
891
892 residues = self.residues[start : end + 1]
893
894 if clone:
895 residues = [r.clone() for r in residues]
896
897 chain = Chain(self.id, accession=self.accession, name=self.name,
898 type=self.type, residues=residues, molecule_id=self.molecule_id)
899 if chain.secondary_structure:
900 chain.secondary_structure = self.secondary_structure.subregion(start, end)
901 chain._torsion_computed = self._torsion_computed
902
903 return chain
904
905 - def find(self, sequence_number, insertion_code=None):
906 """
907 Get a residue by its original Residue Sequence Number and Insertion Code.
908
909 @param sequence_number: PDB sequence number of the residue
910 @type sequence_number: str
911 @param insertion_code: PDB insertion code of the residue (if any)
912 @type insertion_code: str
913
914 @return: the residue object with such an ID
915 @rtype: L{Residue}
916
917 @raise csb.core.ItemNotFoundError: if no residue with that ID exists
918 """
919 res_id = str(sequence_number).strip()
920
921 if insertion_code is not None:
922 insertion_code = str(insertion_code).strip()
923 res_id += insertion_code
924
925 return self.residues._get_residue(res_id)
926
949
951 """
952 Find the optimal fit between C{self} and C{other}. Return L{SuperimposeInfo}
953 (RotationMatrix, translation Vector and RMSD), such that:
954
955 >>> other.transform(rotation_matrix, translation_vector)
956
957 will result in C{other}'s coordinates superimposed over C{self}.
958
959 @param other: the subject (movable) chain
960 @type other: L{Chain}
961 @param what: a list of atom kinds, e.g. ['CA']
962 @type what: list
963 @param how: fitting method (global or local) - a member of the L{AlignmentTypes} enum
964 @type how: L{csb.core.EnumItem}
965
966 @return: superimposition info object, containing rotation matrix, translation
967 vector and computed RMSD
968 @rtype: L{SuperimposeInfo}
969
970 @raise AlignmentArgumentLengthError: when the lengths of the argument chains differ
971 """
972 if self.length != other.length or self.length < 1:
973 raise AlignmentArgumentLengthError('Both chains must be of the same and positive length')
974
975 x = self.get_coordinates(what)
976 y = other.get_coordinates(what)
977 assert len(x) == len(y)
978
979 if how == AlignmentTypes.Global:
980 r, t = csb.bio.utils.fit(x, y)
981 else:
982 r, t = csb.bio.utils.fit_wellordered(x, y)
983
984 rmsd = csb.bio.utils.rmsd(x, y)
985
986 return SuperimposeInfo(r, t, rmsd=rmsd)
987
989 """
990 Align C{other}'s alpha carbons over self in space and return L{SuperimposeInfo}.
991 Coordinates of C{other} are overwritten in place using the rotation matrix
992 and translation vector in L{SuperimposeInfo}. Alias for::
993
994 R, t = self.superimpose(other, what=['CA'])
995 other.transform(R, t)
996
997 @param other: the subject (movable) chain
998 @type other: L{Chain}
999 @param what: a list of atom kinds, e.g. ['CA']
1000 @type what: list
1001 @param how: fitting method (global or local) - a member of the L{AlignmentTypes} enum
1002 @type how: L{csb.core.EnumItem}
1003
1004 @return: superimposition info object, containing rotation matrix, translation
1005 vector and computed RMSD
1006 @rtype: L{SuperimposeInfo}
1007 """
1008 result = self.superimpose(other, what=what, how=how)
1009 other.transform(result.rotation, result.translation)
1010
1011 return result
1012
1013 - def rmsd(self, other, what=['CA']):
1014 """
1015 Compute the C-alpha RMSD against another chain (assuming equal length).
1016 Chains are superimposed with Least Squares Fit / Singular Value Decomposition.
1017
1018 @param other: the subject (movable) chain
1019 @type other: L{Chain}
1020 @param what: a list of atom kinds, e.g. ['CA']
1021 @type what: list
1022
1023 @return: computed RMSD over the specified atom kinds
1024 @rtype: float
1025 """
1026
1027 if self.length != other.length or self.length < 1:
1028 raise ValueError('Both chains must be of the same and positive length '
1029 '(got {0} and {1})'.format(self.length, other.length))
1030
1031 x = self.get_coordinates(what)
1032 y = other.get_coordinates(what)
1033 assert len(x) == len(y)
1034
1035 return csb.bio.utils.rmsd(x, y)
1036
1038 """
1039 Find the optimal fit between C{self} and C{other}. Return L{SuperimposeInfo}
1040 (RotationMatrix, translation Vector and TM-score), such that:
1041
1042 >>> other.transform(rotation_matrix, translation_vector)
1043
1044 will result in C{other}'s coordinates superimposed over C{self}.
1045
1046 @param other: the subject (movable) chain
1047 @type other: L{Chain}
1048 @param what: a list of atom kinds, e.g. ['CA']
1049 @type what: list
1050 @param how: fitting method (global or local) - a member of the L{AlignmentTypes} enum
1051 @type how: L{csb.core.EnumItem}
1052
1053 @return: superimposition info object, containing rotation matrix, translation
1054 vector and computed TM-score
1055 @rtype: L{SuperimposeInfo}
1056
1057 @raise AlignmentArgumentLengthError: when the lengths of the argument chains differ
1058 """
1059
1060 if self.length != other.length or self.length < 1:
1061 raise ValueError('Both chains must be of the same and positive length')
1062
1063 x = self.get_coordinates(what)
1064 y = other.get_coordinates(what)
1065 assert len(x) == len(y)
1066
1067 L_ini_min = 0
1068 if how == AlignmentTypes.Global:
1069 fit = csb.bio.utils.fit
1070 elif how == AlignmentTypes.Local:
1071 fit = csb.bio.utils.fit_wellordered
1072 else:
1073
1074 fit = csb.bio.utils.fit
1075 L_ini_min = 4
1076
1077 r, t, tm = csb.bio.utils.tm_superimpose(x, y, fit, None, None, L_ini_min)
1078
1079 return SuperimposeInfo(r,t, tm_score=tm)
1080
1081 - def tm_score(self, other, what=['CA']):
1082 """
1083 Compute the C-alpha TM-Score against another chain (assuming equal chain length
1084 and optimal configuration - no fitting is done).
1085
1086 @param other: the subject (movable) chain
1087 @type other: L{Chain}
1088 @param what: a list of atom kinds, e.g. ['CA']
1089 @type what: list
1090
1091 @return: computed TM-Score over the specified atom kinds
1092 @rtype: float
1093 """
1094
1095 if self.length != other.length or self.length < 1:
1096 raise ValueError('Both chains must be of the same and positive length')
1097
1098 x = self.get_coordinates(what)
1099 y = other.get_coordinates(what)
1100 assert len(x) == len(y)
1101
1102 return csb.bio.utils.tm_score(x, y)
1103
1105
1114
1116 if len(self) > 0:
1117 return "<ChainResidues: {0} ... {1}>".format(self[self.start_index], self[self.last_index])
1118 else:
1119 return "<ChainResidues: empty>"
1120
1121 @property
1124
1126 """
1127 Append a new residue to the chain.
1128
1129 @param residue: the new residue
1130 @type residue: L{Residue}
1131
1132 @raise DuplicateResidueIDError: if a residue with the same ID already exists
1133 """
1134 if residue.id and residue.id in self.__lookup:
1135 raise DuplicateResidueIDError('A residue with ID {0} is already defined within the chain'.format(residue.id))
1136 index = super(ChainResiduesCollection, self).append(residue)
1137 residue._container = self
1138 self.__container._torsion_computed = False
1139 self._add(residue)
1140 return index
1141
1143 return id in self.__lookup
1144
1146 if id in self.__lookup:
1147 del self.__lookup[id]
1148
1149 - def _add(self, residue):
1151
1157
1158 -class Residue(csb.core.AbstractNIContainer, AbstractEntity):
1159 """
1160 Base class representing a single residue. Provides a dictionary-like
1161 access to the atoms contained in the residue:
1162
1163 >>> residue['CA']
1164 <Atom [3048]: CA>
1165 >>> residue.atoms['CA']
1166 <Atom [3048]: CA>
1167 >>> residue.items
1168 <iterator of Atom-s>
1169
1170 @param rank: rank of the residue with respect to the chain
1171 @type rank: int
1172 @param type: residue type - a member of any L{SequenceAlphabets}
1173 @type type: L{csb.core.EnumItem}
1174 @param sequence_number: PDB sequence number of the residue
1175 @type sequence_number: str
1176 @param insertion_code: PDB insertion code, if any
1177 @type insertion_code: str
1178 """
1179 - def __init__(self, rank, type, sequence_number=None, insertion_code=None):
1180
1181 self._type = None
1182 self._pdb_name = None
1183 self._rank = int(rank)
1184 self._atoms = ResidueAtomsTable(self)
1185 self._secondary_structure = None
1186 self._torsion = None
1187 self._sequence_number = None
1188 self._insertion_code = None
1189 self._container = None
1190
1191 self.type = type
1192 self.id = sequence_number, insertion_code
1193 self._pdb_name = repr(type)
1194
1195 @property
1198
1200 return '<{1} [{0.rank}]: {0.type!r} {0.id}>'.format(self, self.__class__.__name__)
1201
1202 @property
1204 """
1205 Residue type - a member of any sequence alphabet
1206 @rtype: enum item
1207 """
1208 return self._type
1209 @type.setter
1210 - def type(self, type):
1214
1215 @property
1217 """
1218 Residue's position in the sequence (1-based)
1219 @rtype: int
1220 """
1221 return self._rank
1222
1223 @property
1225 """
1226 Secondary structure element this residue is part of
1227 @rtype: L{SecondaryStructureElement}
1228 """
1229 return self._secondary_structure
1230 @secondary_structure.setter
1235
1236 @property
1238 """
1239 Torsion angles
1240 @rtype: L{TorsionAngles}
1241 """
1242 return self._torsion
1243 @torsion.setter
1248
1249 @property
1251 """
1252 Access residue's atoms by atom name
1253 @rtype: L{ResidueAtomsTable}
1254 """
1255 return self._atoms
1256
1257 @property
1259 for atom in self._atoms:
1260 yield self._atoms[atom]
1261
1262 @property
1264 """
1265 PDB sequence number (if residue.has_structure is True)
1266 @rtype: int
1267 """
1268 return self._sequence_number
1269
1270 @property
1272 """
1273 PDB insertion code (if defined)
1274 @rtype: str
1275 """
1276 return self._insertion_code
1277
1278 @property
1280 """
1281 PDB sequence number [+ insertion code]
1282 @rtype: str
1283 """
1284 if self._sequence_number is None:
1285 return None
1286 elif self._insertion_code is not None:
1287 return str(self._sequence_number) + self._insertion_code
1288 else:
1289 return str(self._sequence_number)
1290 @id.setter
1291 - def id(self, value):
1312
1313 @property
1315 """
1316 True if this residue has any atoms
1317 @rtype: bool
1318 """
1319 return len(self.atoms) > 0
1320
1339
1341
1342 container = self._container
1343 self._container = None
1344 clone = copy.deepcopy(self)
1345 self._container = container
1346
1347 return clone
1348
1349 @staticmethod
1350 - def create(sequence_type, *a, **k):
1351 """
1352 Residue factory method, which returns the proper L{Residue} instance based on
1353 the specified C{sequence_type}. All additional arguments are used to initialize
1354 the subclass by passing them automatically to the underlying constructor.
1355
1356 @param sequence_type: create a Residue of that SequenceType
1357 @type sequence_type: L{csb.core.EnumItem}
1358
1359 @return: a new residue of the proper subclass
1360 @rtype: L{Residue} subclass
1361
1362 @raise ValueError: if the sequence type is not known
1363 """
1364 if sequence_type == SequenceTypes.Protein:
1365 return ProteinResidue(*a, **k)
1366 elif sequence_type == SequenceTypes.NucleicAcid:
1367 return NucleicResidue(*a, **k)
1368 elif sequence_type == SequenceTypes.Unknown:
1369 return UnknownResidue(*a, **k)
1370 else:
1371 raise ValueError(sequence_type)
1372
1374 """
1375 Represents a single amino acid residue.
1376
1377 @param rank: rank of the residue with respect to the chain
1378 @type rank: int
1379 @param type: residue type - a member of
1380 L{csb.bio.sequence.SequenceAlphabets.Protein}
1381 @type type: L{csb.core.EnumItem}
1382 @param sequence_number: PDB sequence number of the residue
1383 @type sequence_number: str
1384 @param insertion_code: PDB insertion code, if any
1385 @type insertion_code: str
1386 """
1387
1388 - def __init__(self, rank, type, sequence_number=None, insertion_code=None):
1402
1404 """
1405 Compute the torsion angles of the current residue with neighboring residues
1406 C{prev_residue} and C{next_residue}.
1407
1408 @param prev_residue: the previous residue in the chain
1409 @type prev_residue: L{Residue}
1410 @param next_residue: the next residue in the chain
1411 @type next_residue: L{Residue}
1412 @param strict: if True, L{Broken3DStructureError} is raised if either C{prev_residue}
1413 or C{next_residue} has a broken structure, else the error is silently
1414 ignored and an empty L{TorsionAngles} instance is created
1415 @type strict: bool
1416
1417 @return: a L{TorsionAngles} object, holding the phi, psi and omega values
1418 @rtype: L{TorsionAngles}
1419
1420 @raise Broken3DStructureError: when a specific atom cannot be found
1421 """
1422 if prev_residue is None and next_residue is None:
1423 raise ValueError('At least one neighboring residue is required to compute the torsion.')
1424
1425 angles = TorsionAngles(None, None, None, units=AngleUnits.Degrees)
1426
1427 for residue in (self, prev_residue, next_residue):
1428 if residue is not None and not residue.has_structure:
1429 if strict:
1430 raise Missing3DStructureError(repr(residue))
1431 elif residue is self:
1432 return angles
1433
1434 try:
1435 n = self._atoms['N'].vector
1436 ca = self._atoms['CA'].vector
1437 c = self._atoms['C'].vector
1438 except csb.core.ItemNotFoundError as missing_atom:
1439 if strict:
1440 raise Broken3DStructureError('Could not retrieve {0} atom from the current residue {1!r}.'.format(
1441 missing_atom, self))
1442 else:
1443 return angles
1444
1445 try:
1446 if prev_residue is not None and prev_residue.has_structure:
1447 prev_c = prev_residue._atoms['C'].vector
1448 angles.phi = csb.numeric.dihedral_angle(prev_c, n, ca, c)
1449 except csb.core.ItemNotFoundError as missing_prevatom:
1450 if strict:
1451 raise Broken3DStructureError('Could not retrieve {0} atom from the i-1 residue {1!r}.'.format(
1452 missing_prevatom, prev_residue))
1453 try:
1454 if next_residue is not None and next_residue.has_structure:
1455 next_n = next_residue._atoms['N'].vector
1456 angles.psi = csb.numeric.dihedral_angle(n, ca, c, next_n)
1457 next_ca = next_residue._atoms['CA'].vector
1458 angles.omega = csb.numeric.dihedral_angle(ca, c, next_n, next_ca)
1459 except csb.core.ItemNotFoundError as missing_nextatom:
1460 if strict:
1461 raise Broken3DStructureError('Could not retrieve {0} atom from the i+1 residue {1!r}.'.format(
1462 missing_nextatom, next_residue))
1463
1464 return angles
1465
1467 """
1468 Represents a single nucleotide residue.
1469
1470 @param rank: rank of the residue with respect to the chain
1471 @type rank: int
1472 @param type: residue type - a member of
1473 L{csb.bio.sequence.SequenceAlphabets.Nucleic}
1474 @type type: L{csb.core.EnumItem}
1475 @param sequence_number: PDB sequence number of the residue
1476 @type sequence_number: str
1477 @param insertion_code: PDB insertion code, if any
1478 @type insertion_code: str
1479 """
1480
1481 - def __init__(self, rank, type, sequence_number=None, insertion_code=None):
1496
1498
1499 - def __init__(self, rank, type, sequence_number=None, insertion_code=None):
1503
1505 """
1506 Represents a collection of atoms. Provides dictionary-like access,
1507 where PDB atom names are used for lookup.
1508 """
1509 - def __init__(self, residue, atoms=None):
1517
1519 if len(self) > 0:
1520 return "<ResidueAtoms: {0}>".format(', '.join(self.keys()))
1521 else:
1522 return "<ResidueAtoms: empty>"
1523
1524 @property
1527
1529 """
1530 Append a new Atom to the catalog.
1531
1532 If the atom has an alternate position, a disordered proxy will be created instead and the
1533 atom will be appended to the L{DisorderedAtom}'s list of children. If a disordered atom
1534 with that name already exists, the atom will be appended to its children only.
1535 If an atom with the same name exists, but it was erroneously not marked as disordered,
1536 that terrible condition will be fixed too.
1537
1538 @param atom: the new atom to append
1539 @type atom: L{Atom}
1540
1541 @raise DuplicateAtomIDError: if an atom with the same sequence number and
1542 insertion code already exists in that residue
1543 """
1544 if atom.residue and atom.residue is not self.__residue:
1545 raise InvalidOperation('This atom is part of another residue')
1546 if atom.alternate or (atom.name in self and isinstance(self[atom.name], DisorderedAtom)):
1547 if atom.name not in self:
1548 atom._residue = self.__residue
1549 dis_atom = DisorderedAtom(atom)
1550 super(ResidueAtomsTable, self).append(dis_atom.name, dis_atom)
1551 else:
1552 if not isinstance(self[atom.name], DisorderedAtom):
1553 buggy_atom = self[atom.name]
1554 assert buggy_atom.alternate in (None, False)
1555 buggy_atom.alternate = True
1556 self.update(atom.name, DisorderedAtom(buggy_atom))
1557 if not atom.alternate:
1558 atom.alternate = True
1559 atom._residue = self.__residue
1560 self[atom.name].append(atom)
1561 else:
1562 if atom.name in self:
1563 raise DuplicateAtomIDError('Atom {0} is already defined for {1}'.format(
1564 atom.name, self.__residue))
1565 else:
1566 super(ResidueAtomsTable, self).append(atom.name, atom)
1567 atom._residue = self.__residue
1568
1569 - def update(self, atom_name, atom):
1570 """
1571 Update the atom with the specified name.
1572
1573 @param atom_name: update key
1574 @type atom_name: str
1575 @param atom: new value for this key
1576 @type atom: L{Atom}
1577
1578 @raise ValueError: if C{atom} has a different name than C{atom_name}
1579 """
1580 if atom.name != atom_name:
1581 raise ValueError("Atom's name differs from the specified key.")
1582 if atom.residue is not self.__residue:
1583 atom._residue = self.__residue
1584
1585 super(ResidueAtomsTable, self)._update({atom_name: atom})
1586
1587 -class Atom(AbstractEntity):
1588 """
1589 Represents a single atom in space.
1590
1591 @param serial_number: atom's UID
1592 @type serial_number: int
1593 @param name: atom's name
1594 @type name: str
1595 @param element: corresponding L{ChemElements}
1596 @type element: L{csb.core.EnumItem}
1597 @param vector: atom's coordinates
1598 @type vector: numpy array
1599 @param alternate: if True, means that this is a wobbling atom with multiple alternative
1600 locations
1601 @type alternate: bool
1602 """
1603 - def __init__(self, serial_number, name, element, vector, alternate=False):
1635
1637 return "<Atom [{0.serial_number}]: {0.name}>".format(self)
1638
1641
1646
1658
1667
1668 @property
1670 """
1671 PDB serial number
1672 @rtype: int
1673 """
1674 return self._serial_number
1675 @serial_number.setter
1680
1681 @property
1683 """
1684 PDB atom name (trimmed)
1685 @rtype: str
1686 """
1687 return self._name
1688
1689 @property
1691 """
1692 Chemical element - a member of L{ChemElements}
1693 @rtype: enum item
1694 """
1695 return self._element
1696
1697 @property
1699 """
1700 Residue instance that owns this atom (if available)
1701 @rtype: L{Residue}
1702 """
1703 return self._residue
1704 @residue.setter
1711
1712 @property
1714 """
1715 Atom's 3D coordinates (x, y, z)
1716 @rtype: numpy.array
1717 """
1718 return self._vector
1719 @vector.setter
1721 if numpy.shape(vector) != (3,):
1722 raise ValueError("Three dimensional vector expected")
1723 self._vector = numpy.array(vector)
1724
1725 @property
1727 """
1728 Alternative location flag
1729 @rtype: str
1730 """
1731 return self._alternate
1732 @alternate.setter
1734 self._alternate = value
1735
1736 @property
1738 """
1739 Temperature factor
1740 @rtype: float
1741 """
1742 return self._bfactor
1743 @bfactor.setter
1745 self._bfactor = value
1746
1747 @property
1749 """
1750 Occupancy number
1751 @rtype: float
1752 """
1753 return self._occupancy
1754 @occupancy.setter
1756 self._occupancy = value
1757
1758 @property
1760 """
1761 Charge
1762 @rtype: int
1763 """
1764 return self._charge
1765 @charge.setter
1767 self._charge = value
1768
1769 @property
1772
1774 """
1775 A wobbling atom, which has alternative locations. Each alternative is represented
1776 as a 'normal' L{Atom}. The atom with a highest occupancy is selected as a representative,
1777 hence a DisorderedAtom behaves as a regular L{Atom} (proxy of the representative) as well
1778 as a collection of Atoms.
1779
1780 @param atom: the first atom to be appended to the collection of alternatives. It
1781 is automatically defined as a representative, until a new atom with
1782 higher occupancy is appended to the collection
1783 @type atom: L{Atom}
1784 """
1785
1794
1796 try:
1797 return object.__getattribute__(self, name)
1798 except AttributeError:
1799 subject = object.__getattribute__(self, '_DisorderedAtom__rep')
1800 return getattr(subject, name)
1801
1803 """
1804 Append a new atom to the collection of alternatives.
1805
1806 @param atom: the new alternative
1807 @type atom: L{Atom}
1808 """
1809 self.__update_rep(atom)
1810 self.__alt[atom.alternate] = atom
1811
1812 super(DisorderedAtom, self).append(atom)
1813
1814 - def find(self, altloc):
1815 """
1816 Retrieve a specific atom by its altloc identifier.
1817
1818 @param altloc: alternative location identifier
1819 @type altloc: str
1820
1821 @rtype: L{Atom}
1822 """
1823 if altloc in self.__alt:
1824 return self.__alt[altloc]
1825 else:
1826 for atom in self:
1827 if atom.alternate == altloc:
1828 return Atom
1829
1830 raise EntityNotFoundError(altloc)
1831
1836
1843
1845 return "<DisorderedAtom: {0.length} alternative locations>".format(self)
1846
1848 """
1849 Describes a structural alignment result.
1850
1851 @type rotation: Numpy Array
1852 @type translation: L{Vector}
1853 @type rmsd: float
1854 """
1855 - def __init__(self, rotation, translation, rmsd=None, tm_score=None):
1856
1857 self.rotation = rotation
1858 self.translation = translation
1859 self.rmsd = rmsd
1860 self.tm_score = tm_score
1861
1863 """
1864 Describes a Secondary Structure Element.
1865
1866 @param start: start position with reference to the chain
1867 @type start: float
1868 @param end: end position with reference to the chain
1869 @type end: float
1870 @param type: element type - a member of the L{SecStructures} enum
1871 @type type: csb.core.EnumItem
1872 @param score: secondary structure prediction confidence, if available
1873 @type score: int
1874
1875 @raise IndexError: if start/end positions are out of range
1876 """
1877 - def __init__(self, start, end, type, score=None):
1878
1879 if not (0 < start <= end):
1880 raise IndexError('Element coordinates are out of range: 1 <= start <= end.')
1881
1882 self._start = None
1883 self._end = None
1884 self._type = None
1885 self._score = None
1886
1887 self.start = start
1888 self.end = end
1889 self.type = type
1890
1891 if score is not None:
1892 self.score = score
1893
1896
1901
1904
1906 return "<{0.type!r}: {0.start}-{0.end}>".format(self)
1907
1908 @property
1910 """
1911 Start position (1-based)
1912 @rtype: int
1913 """
1914 return self._start
1915 @start.setter
1916 - def start(self, value):
1922
1923 @property
1925 """
1926 End position (1-based)
1927 @rtype: int
1928 """
1929 return self._end
1930 @end.setter
1931 - def end(self, value):
1937
1938 @property
1940 """
1941 Secondary structure type - a member of L{SecStructures}
1942 @rtype: enum item
1943 """
1944 return self._type
1945 @type.setter
1946 - def type(self, value):
1952
1953 @property
1955 """
1956 Number of residues covered by this element
1957 @rtype: int
1958 """
1959 return self.end - self.start + 1
1960
1961 @property
1963 """
1964 Secondary structure confidence values for each residue
1965 @rtype: L{CollectionContainer}
1966 """
1967 return self._score
1968 @score.setter
1969 - def score(self, scores):
1974
1976 """
1977 Return True if C{self} overlaps with C{other}.
1978
1979 @type other: L{SecondaryStructureElement}
1980 @rtype: bool
1981 """
1982 this = set(range(self.start, self.end + 1))
1983 that = set(range(other.start, other.end + 1))
1984 return not this.isdisjoint(that)
1985
1986 - def merge(self, other):
1987 """
1988 Merge C{self} and C{other}.
1989
1990 @type other: L{SecondaryStructureElement}
1991
1992 @return: a new secondary structure element
1993 @rtype: L{SecondaryStructureElement}
1994
1995 @bug: confidence scores are lost
1996 """
1997 if not self.overlaps(other):
1998 raise ValueError("Can't merge non-overlapping secondary structures")
1999 elif self.type != other.type:
2000 raise ValueError("Can't merge secondary structures of different type")
2001
2002 start = min(self.start, other.start)
2003 end = max(self.end, other.end)
2004 assert self.type == other.type
2005
2006 return SecondaryStructureElement(start, end, self.type)
2007
2009 """
2010 Dump the element as a string.
2011
2012 @return: string representation of the element
2013 @rtype: str
2014 """
2015 return str(self.type) * self.length
2016
2031
2033 """
2034 Describes the secondary structure of a chain.
2035 Provides an index-based access to the secondary structure elements of the chain.
2036
2037 @param string: a secondary structure string (e.g. a PSI-PRED output)
2038 @type string: str
2039 @param conf_string: secondary structure prediction confidence values, if available
2040 @type conf_string: str
2041 """
2042 - def __init__(self, string=None, conf_string=None):
2052
2055
2068
2069 @staticmethod
2070 - def parse(string, conf_string=None):
2071 """
2072 Parse secondary structure from DSSP/PSI-PRED output string.
2073
2074 @param string: a secondary structure string (e.g. a PSI-PRED output)
2075 @type string: str
2076 @param conf_string: secondary structure prediction confidence values, if available
2077 @type conf_string: str
2078
2079 @return: a list of L{SecondaryStructureElement}s.
2080 @rtype: list
2081
2082 @raise ValueError: if the confidence string is not of the same length
2083 """
2084 if not isinstance(string, csb.core.string):
2085 raise TypeError(string)
2086
2087 string = ''.join(re.split('\s+', string))
2088 if conf_string is not None:
2089 conf_string = ''.join(re.split('\s+', conf_string))
2090 if not len(string) == len(conf_string):
2091 raise ValueError('The confidence string has unexpected length.')
2092 motifs = [ ]
2093
2094 if not len(string) > 0:
2095 raise ValueError('Empty Secondary Structure string')
2096
2097 currel = string[0]
2098 start = 0
2099
2100 for i, char in enumerate(string + '.'):
2101
2102 if currel != char:
2103 try:
2104 type = csb.core.Enum.parse(SecStructures, currel)
2105 except csb.core.EnumValueError:
2106 raise UnknownSecStructureError(currel)
2107 confidence = None
2108 if conf_string is not None:
2109 confidence = list(conf_string[start : i])
2110 confidence = list(map(int, confidence))
2111 motif = SecondaryStructureElement(start + 1, i, type, confidence)
2112 motifs.append(motif)
2113
2114 currel = char
2115 start = i
2116
2117 return motifs
2118
2119 @property
2121 """
2122 Start position of the leftmost element
2123 @rtype: int
2124 """
2125 return self._minstart
2126
2127 @property
2129 """
2130 End position of the rightmost element
2131 @rtype: int
2132 """
2133 return self._maxend
2134
2136 """
2137 @return: a deep copy of the object
2138 """
2139 return copy.deepcopy(self)
2140
2142 """
2143 Convert to three-state secondary structure (Helix, Strand, Coil).
2144 """
2145 for e in self:
2146 e.simplify()
2147
2149 """
2150 Get back the string representation of the secondary structure.
2151
2152 @return: a string of secondary structure elements
2153 @rtype: str
2154
2155 @bug: [CSB 0000003] If conflicting elements are found at a given rank,
2156 this position is represented as a coil.
2157 """
2158 gap = str(SecStructures.Gap)
2159 coil = str(SecStructures.Coil)
2160
2161 if chain_length is None:
2162 chain_length = max(e.end for e in self)
2163
2164 ss = []
2165
2166 for pos in range(1, chain_length + 1):
2167 elements = self.at(pos)
2168 if len(elements) > 0:
2169 if len(set(e.type for e in elements)) > 1:
2170 ss.append(coil)
2171 else:
2172 ss.append(elements[0].to_string())
2173 else:
2174 ss.append(gap)
2175
2176 return ''.join(ss)
2177
2178 - def at(self, rank, type=None):
2179 """
2180 @return: all secondary structure elements covering the specifid position
2181 @rtype: tuple of L{SecondaryStructureElement}s
2182 """
2183 return self.scan(start=rank, end=rank, filter=type, loose=True, cut=True)
2184
2185 - def scan(self, start, end, filter=None, loose=True, cut=True):
2186 """
2187 Get all secondary structure elements within the specified [start, end] region.
2188
2189 @param start: the start position of the region, 1-based, inclusive
2190 @type start: int
2191 @param end: the end position of the region, 1-based, inclusive
2192 @type end: int
2193 @param filter: return only elements of the specified L{SecStructures} kind
2194 @type filter: L{csb.core.EnumItem}
2195 @param loose: grab all fully or partially matching elements within the region.
2196 if False, return only the elements which strictly reside within
2197 the region
2198 @type loose: bool
2199 @param cut: if an element is partially overlapping with the start..end region,
2200 cut its start and/or end to make it fit into the region. If False,
2201 return the elements with their real lengths
2202 @type cut: bool
2203
2204 @return: a list of deep-copied L{SecondaryStructureElement}s, sorted by their
2205 start position
2206 @rtype: tuple of L{SecondaryStructureElement}s
2207 """
2208 matches = [ ]
2209
2210 for m in self:
2211 if filter and m.type != filter:
2212 continue
2213
2214 if loose:
2215 if start <= m.start <= end or start <= m.end <= end or (m.start <= start and m.end >= end):
2216 partmatch = copy.deepcopy(m)
2217 if cut:
2218 if partmatch.start < start:
2219 partmatch.start = start
2220 if partmatch.end > end:
2221 partmatch.end = end
2222 if partmatch.score:
2223 partmatch.score = partmatch.score[start : end + 1]
2224 matches.append(partmatch)
2225 else:
2226 if m.start >= start and m.end <= end:
2227 matches.append(copy.deepcopy(m))
2228
2229 matches.sort()
2230 return tuple(matches)
2231
2232 - def q3(self, reference, relaxed=True):
2233 """
2234 Compute Q3 score.
2235
2236 @param reference: reference secondary structure
2237 @type reference: L{SecondaryStructure}
2238 @param relaxed: if True, treat gaps as coils
2239 @type relaxed: bool
2240
2241 @return: the percentage of C{reference} residues with identical
2242 3-state secondary structure.
2243 @rtype: float
2244 """
2245
2246 this = self.clone()
2247 this.to_three_state()
2248
2249 ref = reference.clone()
2250 ref.to_three_state()
2251
2252 total = 0
2253 identical = 0
2254
2255 def at(ss, rank):
2256 elements = ss.at(rank)
2257 if len(elements) == 0:
2258 return None
2259 elif len(elements) > 1:
2260 raise ValueError('Flat secondary structure expected')
2261 else:
2262 return elements[0]
2263
2264 for rank in range(ref.start, ref.end + 1):
2265 q = at(this, rank)
2266 s = at(ref, rank)
2267
2268 if s:
2269 if relaxed or s.type != SecStructures.Gap:
2270 total += 1
2271 if q:
2272 if q.type == s.type:
2273 identical += 1
2274 elif relaxed:
2275 pair = set([q.type, s.type])
2276 match = set([SecStructures.Gap, SecStructures.Coil])
2277 if pair.issubset(match):
2278 identical += 1
2279
2280 if total == 0:
2281 return 0.0
2282 else:
2283 return identical * 100.0 / total
2284
2286 """
2287 Same as C{ss.scan(...cut=True)}, but also shift the start-end positions
2288 of all motifs and return a L{SecondaryStructure} instance instead of a list.
2289
2290 @param start: start position of the subregion, with reference to the chain
2291 @type start: int
2292 @param end: start position of the subregion, with reference to the chain
2293 @type end: int
2294
2295 @return: a deep-copy sub-fragment of the original L{SecondaryStructure}
2296 @rtype: L{SecondaryStructure}
2297 """
2298 sec_struct = SecondaryStructure()
2299
2300 for motif in self.scan(start, end, loose=True, cut=True):
2301
2302 motif.start = motif.start - start + 1
2303 motif.end = motif.end - start + 1
2304 if motif.score:
2305 motif.score = list(motif.score)
2306 sec_struct.append(motif)
2307
2308 return sec_struct
2309
2311 """
2312 Describes a collection of torsion angles. Provides 1-based list-like access.
2313
2314 @param items: an initialization list of L{TorsionAngles}
2315 @type items: list
2316 """
2317 - def __init__(self, items=None, start=1):
2321
2323 if len(self) > 0:
2324 return "<TorsionAnglesList: {0} ... {1}>".format(self[self.start_index], self[self.last_index])
2325 else:
2326 return "<TorsionAnglesList: empty>"
2327
2328 @property
2330 """
2331 List of all phi angles
2332 @rtype: list
2333 """
2334 return [a.phi for a in self]
2335
2336 @property
2338 """
2339 List of all psi angles
2340 @rtype: list
2341 """
2342 return [a.psi for a in self]
2343
2344 @property
2346 """
2347 List of all omega angles
2348 @rtype: list
2349 """
2350 return [a.omega for a in self]
2351
2354
2355 - def rmsd(self, other):
2356 """
2357 Calculate the Circular RSMD against another TorsionAnglesCollection.
2358
2359 @param other: subject (right-hand-term)
2360 @type other: L{TorsionAnglesCollection}
2361
2362 @return: RMSD based on torsion angles
2363 @rtype: float
2364
2365 @raise Broken3DStructureError: on discontinuous torsion angle collections
2366 (phi and psi values are still allowed to be absent at the termini)
2367 @raise ValueError: on mismatching torsion angles collection lengths
2368 """
2369 if len(self) != len(other) or len(self) < 1:
2370 raise ValueError('Both collections must be of the same and positive length')
2371
2372 length = len(self)
2373 query, subject = [], []
2374
2375 for n, (q, s) in enumerate(zip(self, other), start=1):
2376
2377 q = q.copy()
2378 q.to_radians()
2379
2380 s = s.copy()
2381 s.to_radians()
2382
2383 if q.phi is None or s.phi is None:
2384 if n == 1:
2385 q.phi = s.phi = 0.0
2386 else:
2387 raise Broken3DStructureError('Discontinuous torsion angles collection at {0}'.format(n))
2388
2389 if q.psi is None or s.psi is None:
2390 if n == length:
2391 q.psi = s.psi = 0.0
2392 else:
2393 raise Broken3DStructureError('Discontinuous torsion angles collection at {0}'.format(n))
2394
2395 query.append([q.phi, q.psi])
2396 subject.append([s.phi, s.psi])
2397
2398 return csb.bio.utils.torsion_rmsd(numpy.array(query), numpy.array(subject))
2399
2401 """
2402 Describes a collection of phi, psi and omega backbone torsion angles.
2403
2404 It is assumed that the supplied values are either None, or fitting into
2405 the range of [-180, +180] for AngleUnites.Degrees and [0, 2pi] for Radians.
2406
2407 @param phi: phi angle value in C{units}
2408 @type phi: float
2409 @param psi: psi angle value in C{units}
2410 @type psi: float
2411 @param omega: omega angle value in C{units}
2412 @type omega: float
2413 @param units: any of L{AngleUnits}'s enum members
2414 @type units: L{csb.core.EnumItem}
2415
2416 @raise ValueError: on invalid/unknown units
2417 """
2418
2440
2442 return "<TorsionAngles: phi={0.phi}, psi={0.psi}, omega={0.omega}>".format(self)
2443
2446
2448 return self.phi is not None \
2449 or self.psi is not None \
2450 or self.omega is not None
2451
2452 @property
2454 """
2455 Current torsion angle units - a member of L{AngleUnits}
2456 @rtype: enum item
2457 """
2458 return self._units
2459
2460 @property
2463 @phi.setter
2464 - def phi(self, phi):
2467
2468 @property
2471 @psi.setter
2472 - def psi(self, psi):
2475
2476 @property
2479 @omega.setter
2480 - def omega(self, omega):
2483
2489
2506
2507
2524
2525 @staticmethod
2527 """
2528 Check the value of a torsion angle expressed in the specified units.
2529 """
2530 if angle is None:
2531 return
2532 elif units == AngleUnits.Degrees:
2533 if not (-180 <= angle <= 180):
2534 raise ValueError('Torsion angle {0} is out of range -180..180'.format(angle))
2535 elif units == AngleUnits.Radians:
2536 if not (0 <= angle <= (2 * math.pi)):
2537 raise ValueError('Torsion angle {0} is out of range 0..2pi'.format(angle))
2538 else:
2539 raise ValueError('Unknown angle unit type {0}'.format(units))
2540
2541 @staticmethod
2543 """
2544 Convert a torsion angle value, expressed in degrees, to radians.
2545 Negative angles are converted to their positive counterparts: rad(ang + 360deg).
2546
2547 Return the calculated value in the range of [0, 2pi] radians.
2548 """
2549 TorsionAngles.check_angle(angle, AngleUnits.Degrees)
2550
2551 if angle is not None:
2552 if angle < 0:
2553 angle += 360.
2554 angle = math.radians(angle)
2555 return angle
2556
2557 @staticmethod
2559 """
2560 Convert a torsion angle value, expressed in radians, to degrees.
2561 Negative angles are not accepted, it is assumed that negative torsion angles have been
2562 converted to their ang+2pi counterparts beforehand.
2563
2564 Return the calculated value in the range of [-180, +180] degrees.
2565 """
2566 TorsionAngles.check_angle(angle, AngleUnits.Radians)
2567
2568 if angle is not None:
2569 if angle > math.pi:
2570 angle = -((2. * math.pi) - angle)
2571 angle = math.degrees(angle)
2572
2573 return angle
2574