rev |
line source |
me@261
|
1 import sys |
bnagaev@357
|
2 import re |
me@261
|
3 |
me@315
|
4 import util |
me@284
|
5 import fasta |
me@260
|
6 |
dendik@382
|
7 # import this very module as means of having all related classes in one place |
dendik@382
|
8 import base |
dendik@382
|
9 |
me@306
|
10 default_gaps = set((".", "-", "~")) |
me@306
|
11 """Set of characters to recoginze as gaps when parsing alignment.""" |
me@306
|
12 |
me@328
|
13 class Monomer(object): |
me@328
|
14 """Monomer object.""" |
me@260
|
15 |
me@328
|
16 type = None |
me@328
|
17 """Either of 'dna', 'rna', 'protein'.""" |
me@260
|
18 |
dendik@382
|
19 types = base |
dendik@382
|
20 """Mapping of related types. SHOULD be redefined in subclasses.""" |
dendik@382
|
21 |
me@260
|
22 by_code1 = {} |
me@328
|
23 """A mapping from 1-letter code to Monomer subclass.""" |
me@328
|
24 |
me@260
|
25 by_code3 = {} |
me@328
|
26 """A mapping from 3-letter code to Monomer subclass.""" |
me@328
|
27 |
me@260
|
28 by_name = {} |
me@328
|
29 """A mapping from full monomer name to Monomer subclass.""" |
me@260
|
30 |
me@260
|
31 @classmethod |
me@328
|
32 def _subclass(cls, name='', code1='', code3='', is_modified=False): |
me@328
|
33 """Create new subclass of Monomer for given monomer type.""" |
me@328
|
34 class TheMonomer(cls): |
me@328
|
35 pass |
me@328
|
36 name = name.strip().capitalize() |
me@328
|
37 code1 = code1.upper() |
me@328
|
38 code3 = code3.upper() |
bnagaev@357
|
39 TheMonomer.__name__ = re.sub(r"[^\w]", "_", name) |
me@328
|
40 TheMonomer.name = name |
me@328
|
41 TheMonomer.code1 = code1 |
me@328
|
42 TheMonomer.code3 = code3 |
me@328
|
43 TheMonomer.is_modified = is_modified |
me@328
|
44 if not is_modified: |
me@328
|
45 cls.by_code1[code1] = TheMonomer |
me@328
|
46 cls.by_code3[code3] = TheMonomer |
me@328
|
47 cls.by_name[name] = TheMonomer |
me@328
|
48 # We duplicate distinguished long names into Monomer itself, so that we |
me@328
|
49 # can use Monomer.from_code3 to create the relevant type of monomer. |
me@328
|
50 Monomer.by_code3[code3] = TheMonomer |
me@328
|
51 Monomer.by_name[name] = TheMonomer |
me@260
|
52 |
me@328
|
53 @classmethod |
me@353
|
54 def _initialize(cls, codes=None): |
me@328
|
55 """Create all relevant subclasses of Monomer.""" |
me@328
|
56 # NB. The table uses letters d, r, p for types, |
me@328
|
57 # while we use full words; hence, we compare by first letter |
bnagaev@378
|
58 for code1, is_modified, code3, name in codes: |
bnagaev@378
|
59 cls._subclass(name, code1, code3, is_modified) |
me@260
|
60 |
me@260
|
61 @classmethod |
me@260
|
62 def from_code1(cls, code1): |
me@328
|
63 """Create new monomer from 1-letter code.""" |
me@328
|
64 return cls.by_code1[code1.upper()]() |
me@260
|
65 |
me@260
|
66 @classmethod |
me@260
|
67 def from_code3(cls, code3): |
me@328
|
68 """Create new monomer from 3-letter code.""" |
me@328
|
69 return cls.by_code3[code3.upper()]() |
me@260
|
70 |
me@260
|
71 @classmethod |
me@260
|
72 def from_name(cls, name): |
me@328
|
73 """Create new monomer from full name.""" |
me@328
|
74 return cls.by_name[name.strip().capitalize()]() |
me@260
|
75 |
me@329
|
76 def __repr__(self): |
me@329
|
77 return '<Monomer %s>' % self.code3 |
me@329
|
78 |
me@329
|
79 def __str__(self): |
me@329
|
80 """Returns one-letter code""" |
me@329
|
81 return self.code1 |
me@329
|
82 |
me@260
|
83 def __eq__(self, other): |
me@328
|
84 """Monomers within same monomer type are compared by code1.""" |
me@328
|
85 assert self.type == other.type |
me@328
|
86 return self.code1 == other.code1 |
bnagaev@239
|
87 |
bnagaev@239
|
88 class Sequence(list): |
me@274
|
89 """Sequence of Monomers. |
bnagaev@243
|
90 |
me@274
|
91 This behaves like list of monomer objects. In addition to standard list |
me@274
|
92 behaviour, Sequence has the following attributes: |
me@270
|
93 |
me@274
|
94 * name -- str with the name of the sequence |
me@274
|
95 * description -- str with description of the sequence |
me@274
|
96 * source -- str denoting source of the sequence |
me@266
|
97 |
me@274
|
98 Any of them may be empty (i.e. hold empty string) |
me@274
|
99 """ |
me@270
|
100 |
dendik@382
|
101 types = base |
dendik@382
|
102 """Mapping of related types. SHOULD be redefined in subclasses.""" |
me@270
|
103 |
me@275
|
104 name = '' |
me@275
|
105 description = '' |
me@275
|
106 source = '' |
me@275
|
107 |
me@347
|
108 @classmethod |
me@347
|
109 def from_monomers(cls, monomers=[], name=None, description=None, source=None): |
me@347
|
110 """Create sequence from a list of monomer objecst.""" |
bnagaev@378
|
111 result = cls(monomers) |
me@275
|
112 if name: |
me@347
|
113 result.name = name |
me@275
|
114 if description: |
me@347
|
115 result.description = description |
me@275
|
116 if source: |
me@347
|
117 result.source = source |
me@347
|
118 return result |
me@347
|
119 |
me@347
|
120 @classmethod |
me@347
|
121 def from_string(cls, string, name='', description='', source=''): |
me@347
|
122 """Create sequences from string of one-letter codes.""" |
dendik@382
|
123 monomer = cls.types.Monomer.from_code1 |
me@347
|
124 monomers = [monomer(letter) for letter in string] |
me@347
|
125 return cls.from_monomers(monomers, name, description, source) |
me@270
|
126 |
me@329
|
127 def __repr__(self): |
me@329
|
128 return '<Sequence %s>' % str(self) |
me@329
|
129 |
me@262
|
130 def __str__(self): |
me@329
|
131 """Returns sequence of one-letter codes.""" |
me@275
|
132 return ''.join(monomer.code1 for monomer in self) |
me@270
|
133 |
me@316
|
134 def __hash__(self): |
me@316
|
135 """Hash sequence by identity.""" |
me@316
|
136 return id(self) |
me@316
|
137 |
me@295
|
138 class Alignment(object): |
me@295
|
139 """Alignment. It is a list of Columns.""" |
bnagaev@249
|
140 |
dendik@382
|
141 types = base |
dendik@382
|
142 """Mapping of related types. SHOULD be redefined in subclasses.""" |
me@288
|
143 |
me@289
|
144 sequences = None |
me@289
|
145 """Ordered list of sequences in alignment. Read, but DO NOT FIDDLE!""" |
bnagaev@249
|
146 |
me@287
|
147 def __init__(self): |
me@287
|
148 """Initialize empty alignment.""" |
me@287
|
149 self.sequences = [] |
me@295
|
150 self.columns = [] |
me@282
|
151 |
me@362
|
152 # Alignment grow & IO methods |
me@299
|
153 # ============================== |
me@299
|
154 |
me@294
|
155 def append_sequence(self, sequence): |
me@365
|
156 """Add sequence to alignment. Return self. |
me@294
|
157 |
me@294
|
158 If sequence is too short, pad it with gaps on the right. |
me@294
|
159 """ |
me@294
|
160 self.sequences.append(sequence) |
me@294
|
161 for i, monomer in enumerate(sequence): |
me@366
|
162 self._column_at(i)[sequence] = monomer |
me@365
|
163 return self |
me@294
|
164 |
me@364
|
165 def append_row_from_string(self, string, |
me@364
|
166 name='', description='', source='', gaps=default_gaps): |
me@364
|
167 """Add row from a string of one-letter codes and gaps. Return self.""" |
dendik@382
|
168 Sequence = self.types.Sequence |
me@306
|
169 not_gap = lambda (i, char): char not in gaps |
me@349
|
170 without_gaps = util.remove_each(string, gaps) |
me@321
|
171 sequence = Sequence.from_string(without_gaps, name, description, source) |
me@303
|
172 # The following line has some simple magic: |
me@303
|
173 # 1. attach natural numbers to monomers |
me@303
|
174 # 2. delete gaps |
me@303
|
175 # 3. attach numbers again |
me@303
|
176 # This way we have a pair of numbers attached to monomer: |
me@303
|
177 # - it's position in alignment (the first attached number, j) |
me@303
|
178 # - it's position in sequence (the second attached number, i) |
me@349
|
179 for i, (j, char) in enumerate(filter(not_gap, enumerate(string))): |
me@366
|
180 self._column_at(j)[sequence] = sequence[i] |
me@287
|
181 self.sequences.append(sequence) |
me@364
|
182 return self |
me@287
|
183 |
me@366
|
184 def _column_at(self, n): |
me@366
|
185 """Return column by index. Create new columns if required.""" |
me@302
|
186 for i in range(len(self.columns), n + 1): |
me@302
|
187 self.columns.append(Column()) |
me@302
|
188 return self.columns[n] |
me@302
|
189 |
me@362
|
190 def append_file(self, file, format='fasta', gaps=default_gaps): |
me@365
|
191 """Append sequences from file to alignment. Return self. |
me@299
|
192 |
me@362
|
193 If sequences in file have gaps (detected as characters belonging to |
me@362
|
194 `gaps` set), treat them accordingly. |
me@362
|
195 """ |
me@367
|
196 assert format == 'fasta', "We don't support other formats yet" |
me@313
|
197 for (name, description, body) in fasta.parse_file(file): |
bnagaev@378
|
198 self.append_row_from_string(body, name, description, file.name, gaps) |
me@287
|
199 return self |
bnagaev@249
|
200 |
me@367
|
201 def to_file(self, file, format='fasta'): |
me@292
|
202 """Write alignment in FASTA file as sequences with gaps.""" |
me@367
|
203 assert format == "fasta", "We don't support other formats yet" |
me@292
|
204 def char(monomer): |
me@292
|
205 if monomer: |
me@292
|
206 return monomer.code1 |
me@292
|
207 return "-" |
me@292
|
208 for row in self.rows_as_lists(): |
me@292
|
209 seq = row.sequence |
me@292
|
210 line = "".join(map(char, row)) |
me@292
|
211 fasta.save_file(file, line, seq.name, seq.description) |
me@292
|
212 |
me@299
|
213 # Data access methods for alignment |
me@299
|
214 # ================================= |
me@299
|
215 |
me@299
|
216 def rows(self): |
me@299
|
217 """Return list of rows (temporary objects) in alignment. |
me@299
|
218 |
me@299
|
219 Each row is a dictionary of { column : monomer }. |
me@363
|
220 |
me@299
|
221 For gap positions there is no key for the column in row. |
me@299
|
222 |
me@299
|
223 Each row has attribute `sequence` pointing to the sequence the row is |
me@299
|
224 describing. |
me@299
|
225 |
me@299
|
226 Modifications of row have no effect on the alignment. |
me@299
|
227 """ |
me@299
|
228 # For now, the function returns a list rather than iterator. |
me@299
|
229 # It is yet to see, whether memory performance here becomes critical, |
me@299
|
230 # or is random access useful. |
me@299
|
231 rows = [] |
me@299
|
232 for sequence in self.sequences: |
me@299
|
233 row = util.UserDict() |
me@299
|
234 row.sequence = sequence |
me@299
|
235 for column in self.columns: |
me@299
|
236 if sequence in column: |
me@299
|
237 row[column] = column[sequence] |
me@299
|
238 rows.append(row) |
me@299
|
239 return rows |
me@299
|
240 |
me@299
|
241 def rows_as_lists(self): |
me@299
|
242 """Return list of rows (temporary objects) in alignment. |
me@299
|
243 |
me@299
|
244 Each row here is a list of either monomer or None (for gaps). |
me@299
|
245 |
me@299
|
246 Each row has attribute `sequence` pointing to the sequence of row. |
me@299
|
247 |
me@299
|
248 Modifications of row have no effect on the alignment. |
me@299
|
249 """ |
me@299
|
250 rows = [] |
me@299
|
251 for sequence in self.sequences: |
me@299
|
252 row = util.UserList() |
me@299
|
253 row.sequence = sequence |
me@299
|
254 for column in self.columns: |
me@299
|
255 row.append(column.get(sequence)) |
me@299
|
256 rows.append(row) |
me@299
|
257 return rows |
me@299
|
258 |
me@299
|
259 def columns_as_lists(self): |
me@299
|
260 """Return list of columns (temorary objects) in alignment. |
me@299
|
261 |
me@299
|
262 Each column here is a list of either monomer or None (for gaps). |
me@299
|
263 |
me@299
|
264 Items of column are sorted in the same way as alignment.sequences. |
me@299
|
265 |
me@299
|
266 Modifications of column have no effect on the alignment. |
me@299
|
267 """ |
me@299
|
268 columns = [] |
me@299
|
269 for column in self.columns: |
me@299
|
270 col = [] |
me@299
|
271 for sequence in self.sequences: |
me@299
|
272 col.append(column.get(sequence)) |
me@299
|
273 columns.append(col) |
me@299
|
274 return columns |
me@299
|
275 |
me@368
|
276 # Alignment / Block editing methods |
me@368
|
277 # ================================= |
me@368
|
278 |
me@368
|
279 def _flush_row(self, row, whence='left'): |
me@368
|
280 """Helper for `flush`: flush to one side all monomers in one row.""" |
me@368
|
281 row = filter(None, row) |
me@368
|
282 padding = [None] * len(self.columns) |
me@368
|
283 if whence == 'left': |
me@368
|
284 return row + padding |
me@368
|
285 if whence == 'right': |
me@368
|
286 return padding + row |
me@368
|
287 if whence == 'center': |
me@368
|
288 pad_len = (len(self.columns) - len(row)) // 2 |
me@368
|
289 # vvv fix padding for case when length is odd: better have more |
me@368
|
290 pad_len += len(self.columns) - 2 * pad_len |
me@368
|
291 padding = [None] * pad_len |
me@368
|
292 return padding + row + padding |
me@368
|
293 assert True, "whence must be either 'left' or 'right' or 'center'" |
me@368
|
294 |
me@368
|
295 def flush(self, whence='left'): |
me@368
|
296 """Remove all gaps from alignment and flush results to one side. |
me@368
|
297 |
me@368
|
298 `whence` must be one of 'left', 'right' or 'center' |
me@368
|
299 """ |
me@368
|
300 for row in self.rows_as_lists(): |
me@368
|
301 sequence = row.sequence |
me@368
|
302 row = self._flush_row(row, whence) |
me@368
|
303 for monomer, column in zip(row, self.columns): |
me@368
|
304 if monomer: |
me@368
|
305 column[sequence] = monomer |
me@368
|
306 elif sequence in column: |
me@368
|
307 del column[sequence] |
me@368
|
308 |
me@369
|
309 def remove_gap_columns(self): |
me@369
|
310 """Remove all empty columns.""" |
me@369
|
311 for n, column in reversed(enumerate(self.columns)): |
me@369
|
312 if column == {}: |
me@369
|
313 self.columns[n:n+1] = [] |
me@369
|
314 |
me@371
|
315 def _wipe(self): |
me@371
|
316 """Make all positions gaps (but keep sequences intact).""" |
me@371
|
317 for column in self.columns: |
bnagaev@378
|
318 for sequence in list(column.keys()): |
me@371
|
319 del column[sequence] |
me@371
|
320 |
me@372
|
321 def _merge(self, dst, new, merge): |
me@373
|
322 """Replace contents of `dst` with those of `new`. |
me@372
|
323 |
me@372
|
324 Replace contents of elements using function `merge(dst_el, new_le)`. |
me@372
|
325 """ |
dendik@381
|
326 map(merge, zip(dst, new)) |
me@372
|
327 dst[len(dst):] = new[len(dst):] |
me@372
|
328 del dst[len(new):] |
me@371
|
329 |
me@373
|
330 def _replace_sequence_contents(self, new, copy_descriptions): |
me@373
|
331 """Replace contents of sequences with those of `new` alignment.""" |
me@371
|
332 # XXX: we manually copy sequence contents here |
me@372
|
333 # XXX: we only copy, overlapping parts and link to the rest |
me@372
|
334 def merge_monomers(dst, new): |
me@372
|
335 dst.__class__ = new.__class__ |
me@372
|
336 def merge_sequences(dst, new): |
me@373
|
337 if copy_descriptions: |
me@373
|
338 vars(dst).update(vars(new)) |
me@372
|
339 self._merge(dst, new, merge_monomers) |
me@372
|
340 self._merge(self.sequences, new.sequences, merge_sequences) |
me@371
|
341 |
me@371
|
342 def _replace_column_contents(self, new): |
me@373
|
343 """Replace column contents with those of `new` alignment. |
me@371
|
344 |
me@373
|
345 Synonym: copy gap patterns from `new` to `self`. |
me@372
|
346 |
me@373
|
347 `self.sequences` and `new.sequences` should have the same contents. |
me@371
|
348 """ |
me@371
|
349 self._wipe() |
me@371
|
350 not_gap = lambda (a,b): a != None |
me@371
|
351 for sequence, new_row in zip(self.sequences, new.rows_as_lists()): |
me@371
|
352 assert len(sequence) == len(new_row.sequence) |
me@371
|
353 zipped = zip(sequence, filter(not_gap, enumerate(new_row))) |
me@371
|
354 for monomer, (i, _) in zipped: |
me@371
|
355 self._column_at(i)[sequence] = monomer |
me@371
|
356 |
me@373
|
357 def _replace_contents(self, new, copy_descriptions, copy_contents): |
me@371
|
358 """Replace alignment contents with those of other alignment.""" |
me@373
|
359 if copy_contents: |
me@373
|
360 self._replace_sequence_contents(new, copy_descriptions) |
bnagaev@378
|
361 self._replace_column_contents(new) |
me@371
|
362 |
me@373
|
363 def process(self, function, copy_descriptions=True, copy_contents=True): |
me@371
|
364 """Apply function to the alignment (or block); inject results back. |
me@371
|
365 |
me@373
|
366 - `function(block)` must return block with same line order. |
me@373
|
367 - if `copy_descriptions` is False, ignore new sequence names. |
me@373
|
368 - if `copy_contents` is False, don't copy sequence contents too. |
dendik@380
|
369 |
dendik@380
|
370 `function` (object) may have attributes `copy_descriptions` and |
dendik@380
|
371 `copy_contents`, which override the same named arguments. |
me@371
|
372 """ |
me@371
|
373 new = function(self) |
dendik@380
|
374 if hasattr(function, 'copy_descriptions'): |
dendik@380
|
375 copy_descriptions = function.copy_descriptions |
dendik@380
|
376 if hasattr(function, 'copy_contents'): |
dendik@380
|
377 copy_contents = function.copy_contents |
me@373
|
378 self._replace_contents(new, copy_descriptions, copy_contents) |
me@371
|
379 |
me@300
|
380 class Column(dict): |
me@300
|
381 """Column of alignment. |
me@300
|
382 |
me@300
|
383 Column is a dict of { sequence : monomer }. |
me@300
|
384 |
me@300
|
385 For sequences that have gaps in current row, given key is not present in |
me@300
|
386 the column. |
me@300
|
387 """ |
me@325
|
388 |
dendik@382
|
389 types = base |
dendik@382
|
390 """Mapping of related types. SHOULD be redefined in subclasses.""" |
dendik@382
|
391 |
me@325
|
392 def __hash__(self): |
me@325
|
393 """Return hash by identity.""" |
me@325
|
394 return id(self) |
me@300
|
395 |
me@317
|
396 class Block(Alignment): |
me@307
|
397 """Block of alignment. |
me@301
|
398 |
me@307
|
399 Block is intersection of a set of columns & a set of rows. Most of blocks |
me@307
|
400 look like rectangular part of alignment if you shuffle alignment rows the |
me@307
|
401 right way. |
me@261
|
402 """ |
me@270
|
403 |
me@307
|
404 alignment = None |
me@307
|
405 """Alignment the block belongs to.""" |
me@270
|
406 |
me@307
|
407 sequences = () |
me@307
|
408 """List of sequences in block.""" |
me@307
|
409 |
me@307
|
410 columns = () |
me@307
|
411 """List of columns in block.""" |
me@307
|
412 |
me@317
|
413 @classmethod |
me@317
|
414 def from_alignment(cls, alignment, sequences=None, columns=None): |
me@307
|
415 """Build new block from alignment. |
me@307
|
416 |
me@307
|
417 If sequences are not given, the block uses all sequences in alignment. |
me@307
|
418 |
me@307
|
419 If columns are not given, the block uses all columns in alignment. |
me@307
|
420 |
me@307
|
421 In both cases we use exactly the list used in alignment, thus, if new |
me@307
|
422 sequences or columns are added to alignment, the block tracks this too. |
me@261
|
423 """ |
me@307
|
424 if sequences is None: |
me@307
|
425 sequences = alignment.sequences |
me@318
|
426 if columns is None: |
me@307
|
427 columns = alignment.columns |
me@320
|
428 block = cls() |
me@320
|
429 block.alignment = alignment |
me@320
|
430 block.sequences = sequences |
me@320
|
431 block.columns = columns |
me@320
|
432 return block |
me@270
|
433 |
me@260
|
434 # vim: set ts=4 sts=4 sw=4 et: |