1
2
3
4
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
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
30 rl1=Selection.unfold_entities(m1, 'R')
31 rl2=Selection.unfold_entities(m2, 'R')
32
33 p1=0
34 p2=0
35
36 map12={}
37 map21={}
38
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
46 while 1:
47
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
57 while 1:
58
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
68 map12[r1]=r2
69 if r2:
70
71 map21[r2]=r1
72
73 duos.append((r1, r2))
74 self.map12=map12
75 self.map21=map21
76 self.duos=duos
77
83
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
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
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
119 p=PDBParser()
120 s1=p.get_structure('1', pdb_file1)
121 p=PDBParser()
122 s2=p.get_structure('2', pdb_file2)
123
124
125 m1=s1[0]
126 m2=s2[0]
127
128 al=StructureAlignment(fa, m1, m2)
129
130
131 for (r1,r2) in al.get_iterator():
132 print r1, r2
133