1
2
3
4
5
6 """Map the residues of two structures to each other based on a FASTA alignment
7 file.
8 """
9
10 from Bio.SCOP.Raf import to_one_letter_code
11
12 from Bio.PDB import Selection
13 from Bio.PDB.Polypeptide import is_aa
14
15
17 """
18 This class aligns two structures based on an alignment of their
19 sequences.
20 """
21 - def __init__(self, fasta_align, m1, m2, si=0, sj=1):
22 """
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 l=fasta_align.get_alignment_length()
29 s1=fasta_align.get_seq_by_num(si)
30 s2=fasta_align.get_seq_by_num(sj)
31
32 rl1=Selection.unfold_entities(m1, 'R')
33 rl2=Selection.unfold_entities(m2, 'R')
34
35 p1=0
36 p2=0
37
38 map12={}
39 map21={}
40
41 duos=[]
42 for i in range(0, l):
43 column=fasta_align.get_column(i)
44 aa1=column[si]
45 aa2=column[sj]
46 if aa1!="-":
47
48 while 1:
49
50 r1=rl1[p1]
51 p1=p1+1
52 if is_aa(r1):
53 break
54 self._test_equivalence(r1, aa1)
55 else:
56 r1=None
57 if aa2!="-":
58
59 while 1:
60
61 r2=rl2[p2]
62 p2=p2+1
63 if is_aa(r2):
64 break
65 self._test_equivalence(r2, aa2)
66 else:
67 r2=None
68 if r1:
69
70 map12[r1]=r2
71 if r2:
72
73 map21[r2]=r1
74
75 duos.append((r1, r2))
76 self.map12=map12
77 self.map21=map21
78 self.duos=duos
79
85
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
94 """
95 Iterator over all residue pairs.
96 """
97 for i in range(0, len(self.duos)):
98 yield self.duos[i]
99
100
101 if __name__=="__main__":
102 import sys
103 from Bio.Alphabet import generic_protein
104 from Bio import AlignIO
105 from Bio.PDB import PDBParser
106
107 if len(sys.argv) != 4:
108 print "Expects three arguments,"
109 print " - FASTA alignment filename (expect two sequences)"
110 print " - PDB file one"
111 print " - PDB file two"
112 sys.exit()
113
114
115 fa=AlignIO.read(open(sys.argv[1]), "fasta", generic_protein)
116
117 pdb_file1=sys.argv[2]
118 pdb_file2=sys.argv[3]
119
120
121 p=PDBParser()
122 s1=p.get_structure('1', pdb_file1)
123 p=PDBParser()
124 s2=p.get_structure('2', pdb_file2)
125
126
127 m1=s1[0]
128 m2=s2[0]
129
130 al=StructureAlignment(fa, m1, m2)
131
132
133 for (r1,r2) in al.get_iterator():
134 print r1, r2
135