Package Bio :: Package PDB :: Module StructureAlignment'
[hide private]
[frames] | no frames]

Source Code for Module Bio.PDB.StructureAlignment'

 1  # Copyright (C) 2002, Thomas Hamelryck (thamelry@binf.ku.dk) 
 2  # This code is part of the Biopython distribution and governed by its 
 3  # license.  Please see the LICENSE file that should have been included 
 4  # as part of this package. 
 5   
 6  """Map residues of two structures to each other based on a FASTA alignment.""" 
 7   
 8  from __future__ import print_function 
 9   
10  from Bio.Data import SCOPData 
11   
12  from Bio.PDB import Selection 
13  from Bio.PDB.Polypeptide import is_aa 
14   
15   
16 -class StructureAlignment(object):
17 """Class to align two structures based on an alignment of their sequences.""" 18
19 - def __init__(self, fasta_align, m1, m2, si=0, sj=1):
20 """Initialize. 21 22 Attributes: 23 - fasta_align - Alignment object 24 - m1, m2 - two models 25 - si, sj - the sequences in the Alignment object that 26 correspond to the structures 27 28 """ 29 l = fasta_align.get_alignment_length() 30 # Get the residues in the models 31 rl1 = Selection.unfold_entities(m1, 'R') 32 rl2 = Selection.unfold_entities(m2, 'R') 33 # Residue positions 34 p1 = 0 35 p2 = 0 36 # Map equivalent residues to each other 37 map12 = {} 38 map21 = {} 39 # List of residue pairs (None if -) 40 duos = [] 41 for i in range(0, l): 42 column = fasta_align[:, i] 43 aa1 = column[si] 44 aa2 = column[sj] 45 if aa1 != "-": 46 # Position in seq1 is not - 47 while True: 48 # Loop until an aa is found 49 r1 = rl1[p1] 50 p1 = p1 + 1 51 if is_aa(r1): 52 break 53 self._test_equivalence(r1, aa1) 54 else: 55 r1 = None 56 if aa2 != "-": 57 # Position in seq2 is not - 58 while True: 59 # Loop until an aa is found 60 r2 = rl2[p2] 61 p2 = p2 + 1 62 if is_aa(r2): 63 break 64 self._test_equivalence(r2, aa2) 65 else: 66 r2 = None 67 if r1: 68 # Map residue in seq1 to its equivalent in seq2 69 map12[r1] = r2 70 if r2: 71 # Map residue in seq2 to its equivalent in seq1 72 map21[r2] = r1 73 # Append aligned pair (r is None if gap) 74 duos.append((r1, r2)) 75 self.map12 = map12 76 self.map21 = map21 77 self.duos = duos
78
79 - def _test_equivalence(self, r1, aa1):
80 """Test if aa in sequence fits aa in structure.""" 81 resname = r1.get_resname() 82 resname = SCOPData.protein_letters_3to1[resname] 83 assert(aa1 == resname)
84
85 - def get_maps(self):
86 """Map residues between the structures. 87 88 Return two dictionaries that map a residue in one structure to 89 the equivealent residue in the other structure. 90 """ 91 return self.map12, self.map21
92
93 - def get_iterator(self):
94 """Iterator over all residue pairs.""" 95 for i in range(0, len(self.duos)): 96 yield self.duos[i]
97