allpy: allpy/base.py annotate

allpy

annotate allpy/base.py @ 347:b6c1eaad5ac6

base: introduced Sequence.from_monomers instead of Sequence() (rels #12)

author	Daniil Alexeyevsky <me.dendik@gmail.com>
date	Thu, 20 Jan 2011 22:06:57 +0300
parents	7880ea44e56b
children	b5b7fcc6c9cb

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