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