rev |
line source |
me@261
|
1 import sys |
me@261
|
2 |
me@315
|
3 import util |
me@284
|
4 import fasta |
me@260
|
5 import data.codes |
me@260
|
6 |
me@306
|
7 default_gaps = set((".", "-", "~")) |
me@306
|
8 """Set of characters to recoginze as gaps when parsing alignment.""" |
me@306
|
9 |
me@260
|
10 class MonomerType(object): |
me@260
|
11 """Class of monomer types. |
me@260
|
12 |
me@260
|
13 Each MonomerType object represents a known monomer type, e.g. Valine, |
me@260
|
14 and is referenced to by each instance of monomer in a given sequence. |
me@260
|
15 |
me@260
|
16 - `name`: full name of monomer type |
me@260
|
17 - `code1`: one-letter code |
me@260
|
18 - `code3`: three-letter code |
me@260
|
19 - `is_modified`: either of True or False |
me@260
|
20 |
me@260
|
21 class atributes: |
me@260
|
22 |
me@260
|
23 - `by_code1`: a mapping from one-letter code to MonomerType object |
me@260
|
24 - `by_code3`: a mapping from three-letter code to MonomerType object |
me@260
|
25 - `by_name`: a mapping from monomer name to MonomerType object |
me@260
|
26 - `instance_type`: class of Monomer objects to use when creating new |
me@260
|
27 objects; this must be redefined in descendent classes |
me@260
|
28 |
me@260
|
29 All of the class attributes MUST be redefined when subclassing. |
me@260
|
30 """ |
me@260
|
31 |
me@260
|
32 by_code1 = {} |
me@260
|
33 by_code3 = {} |
me@260
|
34 by_name = {} |
me@260
|
35 instance_type = None |
me@260
|
36 |
me@260
|
37 def __init__(self, name="", code1="", code3="", is_modified=False): |
me@310
|
38 super(MonomerType, self).__init__() |
me@260
|
39 self.name = name.capitalize() |
me@260
|
40 self.code1 = code1.upper() |
me@260
|
41 self.code3 = code3.upper() |
me@260
|
42 self.is_modified = bool(is_modified) |
me@260
|
43 if not is_modified: |
me@260
|
44 self.by_code1[self.code1] = self |
me@260
|
45 self.by_code3[code3] = self |
me@260
|
46 self.by_name[name] = self |
me@260
|
47 # We duplicate distinguished long names into MonomerType itself, |
me@260
|
48 # so that we can use MonomerType.from_code3 to create the relevant |
me@260
|
49 # type of monomer. |
me@260
|
50 MonomerType.by_code3[code3] = self |
me@260
|
51 MonomerType.by_name[name] = self |
me@260
|
52 |
me@260
|
53 @classmethod |
me@260
|
54 def _initialize(cls, type_letter, codes=data.codes.codes): |
me@260
|
55 """Create all relevant instances of MonomerType. |
me@260
|
56 |
me@260
|
57 `type_letter` is either of: |
me@260
|
58 |
me@260
|
59 - 'p' for protein |
me@260
|
60 - 'd' for DNA |
me@260
|
61 - 'r' for RNA |
me@260
|
62 |
me@260
|
63 `codes` is a table of monomer codes |
me@260
|
64 """ |
me@260
|
65 for type, code1, is_modified, code3, name in codes: |
me@260
|
66 if type == type_letter: |
me@260
|
67 cls(name, code1, code3, is_modified) |
me@260
|
68 |
me@260
|
69 @classmethod |
me@260
|
70 def from_code1(cls, code1): |
me@260
|
71 """Return monomer type by one-letter code.""" |
me@260
|
72 return cls.by_code1[code1.upper()] |
me@260
|
73 |
me@260
|
74 @classmethod |
me@260
|
75 def from_code3(cls, code3): |
me@260
|
76 """Return monomer type by three-letter code.""" |
me@260
|
77 return cls.by_code3[code3.upper()] |
me@260
|
78 |
me@260
|
79 @classmethod |
me@260
|
80 def from_name(cls, name): |
me@260
|
81 """Return monomer type by name.""" |
me@260
|
82 return cls.by_name[name.capitalize()] |
me@260
|
83 |
me@260
|
84 def instance(self): |
me@260
|
85 """Create a new monomer of given type.""" |
me@260
|
86 return self.instance_type(self) |
me@260
|
87 |
me@260
|
88 def __eq__(self, other): |
me@260
|
89 if hasattr(other, "type"): |
me@260
|
90 return self is other.type |
me@260
|
91 return self is other |
me@260
|
92 |
me@260
|
93 class Monomer(object): |
me@260
|
94 """Monomer object. |
me@260
|
95 |
me@260
|
96 attributes: |
me@260
|
97 |
me@260
|
98 - `type`: type of monomer (a MonomerType object) |
me@260
|
99 |
me@282
|
100 class attributes: |
me@282
|
101 |
me@282
|
102 - `monomer_type`: either MonomerType or one of it's subclasses, it is used |
me@282
|
103 when creating new monomers. It SHOULD be redefined when subclassing |
me@282
|
104 Monomer. |
me@260
|
105 """ |
me@260
|
106 monomer_type = MonomerType |
me@260
|
107 |
me@260
|
108 def __init__(self, type): |
me@310
|
109 super(Monomer, self).__init__() |
me@260
|
110 self.type = type |
me@260
|
111 |
me@260
|
112 @classmethod |
me@260
|
113 def from_code1(cls, code1): |
me@260
|
114 return cls(cls.monomer_type.by_code1[code1.upper()]) |
me@260
|
115 |
me@260
|
116 @classmethod |
me@260
|
117 def from_code3(cls, code3): |
me@260
|
118 return cls(cls.monomer_type.by_code3[code3.upper()]) |
me@260
|
119 |
me@260
|
120 @classmethod |
me@260
|
121 def from_name(cls, name): |
me@260
|
122 return cls(cls.monomer_type.by_name[name.capitalize()]) |
me@260
|
123 |
me@260
|
124 def __eq__(self, other): |
me@260
|
125 if hasattr(other, "type"): |
me@260
|
126 return self.type is other.type |
me@260
|
127 return self.type is other |
bnagaev@239
|
128 |
bnagaev@239
|
129 class Sequence(list): |
me@274
|
130 """Sequence of Monomers. |
bnagaev@243
|
131 |
me@274
|
132 This behaves like list of monomer objects. In addition to standard list |
me@274
|
133 behaviour, Sequence has the following attributes: |
me@270
|
134 |
me@274
|
135 * name -- str with the name of the sequence |
me@274
|
136 * description -- str with description of the sequence |
me@274
|
137 * source -- str denoting source of the sequence |
me@266
|
138 |
me@274
|
139 Any of them may be empty (i.e. hold empty string) |
me@275
|
140 |
me@275
|
141 Class attributes: |
me@282
|
142 |
me@275
|
143 * monomer_type -- type of monomers in sequence, must be redefined when |
me@275
|
144 subclassing |
me@274
|
145 """ |
me@270
|
146 |
me@275
|
147 monomer_type = Monomer |
me@270
|
148 |
me@275
|
149 name = '' |
me@275
|
150 description = '' |
me@275
|
151 source = '' |
me@275
|
152 |
me@275
|
153 def __init__(self, sequence=[], name=None, description=None, source=None): |
me@275
|
154 super(Sequence, self).__init__(sequence) |
me@275
|
155 if hasattr(sequence, 'name'): |
me@275
|
156 vars(self).update(vars(sequence)) |
me@275
|
157 if name: |
me@275
|
158 self.name = name |
me@275
|
159 if description: |
me@275
|
160 self.description = description |
me@275
|
161 if source: |
me@275
|
162 self.source = source |
me@270
|
163 |
me@262
|
164 def __str__(self): |
me@275
|
165 """Returns sequence in one-letter code.""" |
me@275
|
166 return ''.join(monomer.code1 for monomer in self) |
me@270
|
167 |
me@273
|
168 @classmethod |
me@273
|
169 def from_string(cls, string, name='', description=''): |
me@273
|
170 """Create sequences from string of one-letter codes.""" |
me@273
|
171 monomer = cls.monomer_type.from_code1 |
me@273
|
172 monomers = [monomer(letter) for letter in string] |
me@273
|
173 return cls(monomers, name, description) |
me@262
|
174 |
me@284
|
175 @classmethod |
me@284
|
176 def from_fasta(cls, file): |
me@284
|
177 """Read sequence from FASTA file. |
me@286
|
178 |
me@284
|
179 File must contain exactly one sequence. |
me@284
|
180 """ |
me@313
|
181 sequence, = fasta.parse_file(file) |
me@313
|
182 name, description, body = sequence |
me@313
|
183 return cls(body, name, description, file.name) |
me@284
|
184 |
me@295
|
185 class Alignment(object): |
me@295
|
186 """Alignment. It is a list of Columns.""" |
bnagaev@249
|
187 |
me@287
|
188 sequence_type = Sequence |
me@289
|
189 """Type of sequences in alignment. SHOULD be redefined when subclassing.""" |
me@288
|
190 |
me@289
|
191 sequences = None |
me@289
|
192 """Ordered list of sequences in alignment. Read, but DO NOT FIDDLE!""" |
bnagaev@249
|
193 |
me@287
|
194 def __init__(self): |
me@287
|
195 """Initialize empty alignment.""" |
me@287
|
196 super(Alignment, self).__init__() |
me@287
|
197 self.sequences = [] |
me@295
|
198 self.columns = [] |
me@282
|
199 |
me@299
|
200 # Alignment modification methods |
me@299
|
201 # ============================== |
me@299
|
202 |
me@294
|
203 def append_sequence(self, sequence): |
me@294
|
204 """Add sequence to alignment. |
me@294
|
205 |
me@294
|
206 If sequence is too short, pad it with gaps on the right. |
me@294
|
207 """ |
me@294
|
208 self.sequences.append(sequence) |
me@294
|
209 for i, monomer in enumerate(sequence): |
me@302
|
210 self.column_at(i)[sequence] = monomer |
me@294
|
211 |
me@306
|
212 def append_gapped_line(self, line, name='', description='', source='', |
me@306
|
213 gaps=default_gaps): |
me@287
|
214 """Add row from a line of one-letter codes and gaps.""" |
me@313
|
215 Sequence = self.sequence_type |
me@306
|
216 not_gap = lambda (i, char): char not in gaps |
me@306
|
217 without_gaps = util.remove_each(line, gaps) |
me@306
|
218 sequence = Sequence(without_gaps, name, description, source) |
me@303
|
219 # The following line has some simple magic: |
me@303
|
220 # 1. attach natural numbers to monomers |
me@303
|
221 # 2. delete gaps |
me@303
|
222 # 3. attach numbers again |
me@303
|
223 # This way we have a pair of numbers attached to monomer: |
me@303
|
224 # - it's position in alignment (the first attached number, j) |
me@303
|
225 # - it's position in sequence (the second attached number, i) |
me@287
|
226 for i, (j, char) in enumerate(filter(not_gap, enumerate(line))): |
me@313
|
227 self.column_at(j)[sequence] = sequence[i] |
me@287
|
228 self.sequences.append(sequence) |
me@287
|
229 |
me@302
|
230 def column_at(self, n): |
me@302
|
231 """Return column by index. Create required new columns if required. |
me@302
|
232 |
me@302
|
233 Do NOT use this method, unless you are sure it is what you want. |
me@302
|
234 """ |
me@302
|
235 for i in range(len(self.columns), n + 1): |
me@302
|
236 self.columns.append(Column()) |
me@302
|
237 return self.columns[n] |
me@302
|
238 |
me@299
|
239 # Alignment IO methods |
me@299
|
240 # ==================== |
me@299
|
241 |
me@287
|
242 @classmethod |
me@306
|
243 def from_fasta(cls, file, gaps=default_gaps): |
me@287
|
244 """Create new alignment from FASTA file.""" |
me@287
|
245 self = cls() |
me@313
|
246 for (name, description, body) in fasta.parse_file(file): |
me@306
|
247 self.append_gapped_line(body, name, description, gaps) |
me@287
|
248 return self |
bnagaev@249
|
249 |
me@292
|
250 def to_fasta(self, file): |
me@292
|
251 """Write alignment in FASTA file as sequences with gaps.""" |
me@292
|
252 def char(monomer): |
me@292
|
253 if monomer: |
me@292
|
254 return monomer.code1 |
me@292
|
255 return "-" |
me@292
|
256 for row in self.rows_as_lists(): |
me@292
|
257 seq = row.sequence |
me@292
|
258 line = "".join(map(char, row)) |
me@292
|
259 fasta.save_file(file, line, seq.name, seq.description) |
me@292
|
260 |
me@299
|
261 # Data access methods for alignment |
me@299
|
262 # ================================= |
me@299
|
263 |
me@299
|
264 def rows(self): |
me@299
|
265 """Return list of rows (temporary objects) in alignment. |
me@299
|
266 |
me@299
|
267 Each row is a dictionary of { column : monomer }. |
me@299
|
268 |
me@299
|
269 For gap positions there is no key for the column in row. |
me@299
|
270 |
me@299
|
271 Each row has attribute `sequence` pointing to the sequence the row is |
me@299
|
272 describing. |
me@299
|
273 |
me@299
|
274 Modifications of row have no effect on the alignment. |
me@299
|
275 """ |
me@299
|
276 # For now, the function returns a list rather than iterator. |
me@299
|
277 # It is yet to see, whether memory performance here becomes critical, |
me@299
|
278 # or is random access useful. |
me@299
|
279 rows = [] |
me@299
|
280 for sequence in self.sequences: |
me@299
|
281 row = util.UserDict() |
me@299
|
282 row.sequence = sequence |
me@299
|
283 for column in self.columns: |
me@299
|
284 if sequence in column: |
me@299
|
285 row[column] = column[sequence] |
me@299
|
286 rows.append(row) |
me@299
|
287 return rows |
me@299
|
288 |
me@299
|
289 def rows_as_lists(self): |
me@299
|
290 """Return list of rows (temporary objects) in alignment. |
me@299
|
291 |
me@299
|
292 Each row here is a list of either monomer or None (for gaps). |
me@299
|
293 |
me@299
|
294 Each row has attribute `sequence` pointing to the sequence of row. |
me@299
|
295 |
me@299
|
296 Modifications of row have no effect on the alignment. |
me@299
|
297 """ |
me@299
|
298 rows = [] |
me@299
|
299 for sequence in self.sequences: |
me@299
|
300 row = util.UserList() |
me@299
|
301 row.sequence = sequence |
me@299
|
302 for column in self.columns: |
me@299
|
303 row.append(column.get(sequence)) |
me@299
|
304 rows.append(row) |
me@299
|
305 return rows |
me@299
|
306 |
me@299
|
307 def columns_as_lists(self): |
me@299
|
308 """Return list of columns (temorary objects) in alignment. |
me@299
|
309 |
me@299
|
310 Each column here is a list of either monomer or None (for gaps). |
me@299
|
311 |
me@299
|
312 Items of column are sorted in the same way as alignment.sequences. |
me@299
|
313 |
me@299
|
314 Modifications of column have no effect on the alignment. |
me@299
|
315 """ |
me@299
|
316 columns = [] |
me@299
|
317 for column in self.columns: |
me@299
|
318 col = [] |
me@299
|
319 for sequence in self.sequences: |
me@299
|
320 col.append(column.get(sequence)) |
me@299
|
321 columns.append(col) |
me@299
|
322 return columns |
me@299
|
323 |
me@300
|
324 class Column(dict): |
me@300
|
325 """Column of alignment. |
me@300
|
326 |
me@300
|
327 Column is a dict of { sequence : monomer }. |
me@300
|
328 |
me@300
|
329 For sequences that have gaps in current row, given key is not present in |
me@300
|
330 the column. |
me@300
|
331 """ |
me@300
|
332 pass |
me@300
|
333 |
me@307
|
334 class Block(object): |
me@307
|
335 """Block of alignment. |
me@301
|
336 |
me@307
|
337 Block is intersection of a set of columns & a set of rows. Most of blocks |
me@307
|
338 look like rectangular part of alignment if you shuffle alignment rows the |
me@307
|
339 right way. |
me@261
|
340 """ |
me@270
|
341 |
me@307
|
342 alignment = None |
me@307
|
343 """Alignment the block belongs to.""" |
me@270
|
344 |
me@307
|
345 sequences = () |
me@307
|
346 """List of sequences in block.""" |
me@307
|
347 |
me@307
|
348 columns = () |
me@307
|
349 """List of columns in block.""" |
me@307
|
350 |
me@307
|
351 def __init__(self, alignment, sequences=None, columns=None): |
me@307
|
352 """Build new block from alignment. |
me@307
|
353 |
me@307
|
354 If sequences are not given, the block uses all sequences in alignment. |
me@307
|
355 |
me@307
|
356 If columns are not given, the block uses all columns in alignment. |
me@307
|
357 |
me@307
|
358 In both cases we use exactly the list used in alignment, thus, if new |
me@307
|
359 sequences or columns are added to alignment, the block tracks this too. |
me@261
|
360 """ |
me@310
|
361 super(Block, self).__init__() |
me@307
|
362 if sequences is None: |
me@307
|
363 sequences = alignment.sequences |
me@307
|
364 if colums is None: |
me@307
|
365 columns = alignment.columns |
me@261
|
366 self.alignment = alignment |
me@261
|
367 self.sequences = sequences |
me@261
|
368 self.positions = positions |
me@270
|
369 |
me@312
|
370 def flush_left(self): |
me@312
|
371 """Move all monomers to the left, gaps to the right within block.""" |
me@312
|
372 padding = [None] * len(self.columns) |
me@312
|
373 for row in self.rows_as_lists(): |
me@312
|
374 sequence = row.sequence |
me@312
|
375 row = filter(None, row) + padding |
me@312
|
376 for monomer, column in zip(row, self.columns): |
me@312
|
377 if monomer: |
me@312
|
378 column[sequence] = monomer |
me@312
|
379 elif sequence in column: |
me@312
|
380 del column[sequence] |
me@312
|
381 |
me@312
|
382 |
me@260
|
383 # vim: set ts=4 sts=4 sw=4 et: |