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  from Bio.PDB import * 
  7   
  8  __doc__=""" 
  9  Map the residues of two structures to each other based on  
 10  a FASTA alignment file. 
 11  """ 
 12   
 13   
14 -class StructureAlignment:
15 """ 16 This class aligns two structures based on an alignment of their 17 sequences. 18 """
19 - def __init__(self, fasta_align, m1, m2, si=0, sj=1):
20 """ 21 fasta_align --- Alignment object 22 m1, m2 --- two models 23 si, sj --- the sequences in the Alignment object that 24 correspond to the structures 25 """ 26 l=fasta_align.get_alignment_length() 27 s1=fasta_align.get_seq_by_num(si) 28 s2=fasta_align.get_seq_by_num(sj) 29 # Get the residues in the models 30 rl1=Selection.unfold_entities(m1, 'R') 31 rl2=Selection.unfold_entities(m2, 'R') 32 # Residue positions 33 p1=0 34 p2=0 35 # Map equivalent residues to each other 36 map12={} 37 map21={} 38 # List of residue pairs (None if -) 39 duos=[] 40 for i in range(0, l): 41 column=fasta_align.get_column(i) 42 aa1=column[si] 43 aa2=column[sj] 44 if aa1!="-": 45 # Position in seq1 is not - 46 while 1: 47 # Loop until an aa is found 48 r1=rl1[p1] 49 p1=p1+1 50 if is_aa(r1): 51 break 52 self._test_equivalence(r1, aa1) 53 else: 54 r1=None 55 if aa2!="-": 56 # Position in seq2 is not - 57 while 1: 58 # Loop until an aa is found 59 r2=rl2[p2] 60 p2=p2+1 61 if is_aa(r2): 62 break 63 self._test_equivalence(r2, aa2) 64 else: 65 r2=None 66 if r1: 67 # Map residue in seq1 to its equivalent in seq2 68 map12[r1]=r2 69 if r2: 70 # Map residue in seq2 to its equivalent in seq1 71 map21[r2]=r1 72 # Append aligned pair (r is None if gap) 73 duos.append((r1, r2)) 74 self.map12=map12 75 self.map21=map21 76 self.duos=duos
77
78 - def _test_equivalence(self, r1, aa1):
79 "Test if aa in sequence fits aa in structure." 80 resname=r1.get_resname() 81 resname=to_one_letter_code[resname] 82 assert(aa1==resname)
83
84 - def get_maps(self):
85 """ 86 Return two dictionaries that map a residue in one structure to 87 the equivealent residue in the other structure. 88 """ 89 return self.map12, self.map21
90
91 - def get_iterator(self):
92 """ 93 Iterator over all residue pairs. 94 """ 95 for i in range(0, len(self.duos)): 96 yield self.duos[i]
97 98 99 if __name__=="__main__": 100 import sys 101 from Bio.Alphabet import generic_protein 102 from Bio import AlignIO 103 from Bio.PDB import * 104 105 if len(sys.argv) != 4 : 106 print "Expects three arguments," 107 print " - FASTA alignment filename (expect two sequences)" 108 print " - PDB file one" 109 print " - PDB file two" 110 sys.exit() 111 112 # The alignment 113 fa=AlignIO.read(open(sys.argv[1]), "fasta", generic_protein) 114 115 pdb_file1=sys.argv[2] 116 pdb_file2=sys.argv[3] 117 118 # The structures 119 p=PDBParser() 120 s1=p.get_structure('1', pdb_file1) 121 p=PDBParser() 122 s2=p.get_structure('2', pdb_file2) 123 124 # Get the models 125 m1=s1[0] 126 m2=s2[0] 127 128 al=StructureAlignment(fa, m1, m2) 129 130 # Print aligned pairs (r is None if gap) 131 for (r1,r2) in al.get_iterator(): 132 print r1, r2 133