allpy: allpy/base.py annotate

allpy

annotate allpy/base.py @ 275:4a2341bc90b1

base.Sequence major cleanup * moved base.Sequence.secstr_has to pdb (it appears that secondary structure it uses is from PDB, which is not a nice thing to do; secondary structure is an entity in itself and must have same names across different sources) * removed stupid base.Sequence.monomers, since base.Sequence is subclass of list * removed all methods that reimplement list behaviour * fixed headers of base.Sequence.from_slice

author	Daniil Alexeyevsky <me.dendik@gmail.com>
date	Wed, 15 Dec 2010 20:06:27 +0300
parents	9ea5099309b8
children	4466baf642df

rev	line source
me@261	1 import sys
me@261	2 import os
me@262	3 import os.path
me@261	4 from tempfile import NamedTemporaryFile
me@262	5 import urllib2
me@261	6
me@261	7 import config
me@261	8 from graph import Graph
me@262	9 from Bio.PDB.DSSP import make_dssp_dict
bnagaev@249	10 from fasta import save_fasta
me@260	11 import data.codes
me@260	12
me@260	13 class MonomerType(object):
me@260	14 """Class of monomer types.
me@260	15
me@260	16 Each MonomerType object represents a known monomer type, e.g. Valine,
me@260	17 and is referenced to by each instance of monomer in a given sequence.
me@260	18
me@260	19 - `name`: full name of monomer type
me@260	20 - `code1`: one-letter code
me@260	21 - `code3`: three-letter code
me@260	22 - `is_modified`: either of True or False
me@260	23
me@260	24 class atributes:
me@260	25
me@260	26 - `by_code1`: a mapping from one-letter code to MonomerType object
me@260	27 - `by_code3`: a mapping from three-letter code to MonomerType object
me@260	28 - `by_name`: a mapping from monomer name to MonomerType object
me@260	29 - `instance_type`: class of Monomer objects to use when creating new
me@260	30 objects; this must be redefined in descendent classes
me@260	31
me@260	32 All of the class attributes MUST be redefined when subclassing.
me@260	33 """
me@260	34
me@260	35 by_code1 = {}
me@260	36 by_code3 = {}
me@260	37 by_name = {}
me@260	38 instance_type = None
me@260	39
me@260	40 def __init__(self, name="", code1="", code3="", is_modified=False):
me@260	41 self.name = name.capitalize()
me@260	42 self.code1 = code1.upper()
me@260	43 self.code3 = code3.upper()
me@260	44 self.is_modified = bool(is_modified)
me@260	45 if not is_modified:
me@260	46 self.by_code1[self.code1] = self
me@260	47 self.by_code3[code3] = self
me@260	48 self.by_name[name] = self
me@260	49 # We duplicate distinguished long names into MonomerType itself,
me@260	50 # so that we can use MonomerType.from_code3 to create the relevant
me@260	51 # type of monomer.
me@260	52 MonomerType.by_code3[code3] = self
me@260	53 MonomerType.by_name[name] = self
me@260	54
me@260	55 @classmethod
me@260	56 def _initialize(cls, type_letter, codes=data.codes.codes):
me@260	57 """Create all relevant instances of MonomerType.
me@260	58
me@260	59 `type_letter` is either of:
me@260	60
me@260	61 - 'p' for protein
me@260	62 - 'd' for DNA
me@260	63 - 'r' for RNA
me@260	64
me@260	65 `codes` is a table of monomer codes
me@260	66 """
me@260	67 for type, code1, is_modified, code3, name in codes:
me@260	68 if type == type_letter:
me@260	69 cls(name, code1, code3, is_modified)
me@260	70
me@260	71 @classmethod
me@260	72 def from_code1(cls, code1):
me@260	73 """Return monomer type by one-letter code."""
me@260	74 return cls.by_code1[code1.upper()]
me@260	75
me@260	76 @classmethod
me@260	77 def from_code3(cls, code3):
me@260	78 """Return monomer type by three-letter code."""
me@260	79 return cls.by_code3[code3.upper()]
me@260	80
me@260	81 @classmethod
me@260	82 def from_name(cls, name):
me@260	83 """Return monomer type by name."""
me@260	84 return cls.by_name[name.capitalize()]
me@260	85
me@260	86 def instance(self):
me@260	87 """Create a new monomer of given type."""
me@260	88 return self.instance_type(self)
me@260	89
me@260	90 def __eq__(self, other):
me@260	91 if hasattr(other, "type"):
me@260	92 return self is other.type
me@260	93 return self is other
me@260	94
me@260	95 class Monomer(object):
me@260	96 """Monomer object.
me@260	97
me@260	98 attributes:
me@260	99
me@260	100 - `type`: type of monomer (a MonomerType object)
me@260	101
me@260	102 class attribute `monomer_type` is MonomerType or either of it's subclasses,
me@260	103 it is used when creating new monomers. It MUST be redefined when subclassing Monomer.
me@260	104 """
me@260	105 monomer_type = MonomerType
me@260	106
me@260	107 def __init__(self, type):
me@260	108 self.type = type
me@260	109
me@260	110 @classmethod
me@260	111 def from_code1(cls, code1):
me@260	112 return cls(cls.monomer_type.by_code1[code1.upper()])
me@260	113
me@260	114 @classmethod
me@260	115 def from_code3(cls, code3):
me@260	116 return cls(cls.monomer_type.by_code3[code3.upper()])
me@260	117
me@260	118 @classmethod
me@260	119 def from_name(cls, name):
me@260	120 return cls(cls.monomer_type.by_name[name.capitalize()])
me@260	121
me@260	122 def __eq__(self, other):
me@260	123 if hasattr(other, "type"):
me@260	124 return self.type is other.type
me@260	125 return self.type is other
bnagaev@239	126
bnagaev@239	127 class Sequence(list):
me@274	128 """Sequence of Monomers.
bnagaev@243	129
me@274	130 This behaves like list of monomer objects. In addition to standard list
me@274	131 behaviour, Sequence has the following attributes:
me@270	132
me@274	133 * name -- str with the name of the sequence
me@274	134 * description -- str with description of the sequence
me@274	135 * source -- str denoting source of the sequence
me@266	136
me@274	137 Any of them may be empty (i.e. hold empty string)
me@275	138
me@275	139 Class attributes:
me@275	140 * monomer_type -- type of monomers in sequence, must be redefined when
me@275	141 subclassing
me@274	142 """
me@270	143
me@275	144 monomer_type = Monomer
me@270	145
me@275	146 name = ''
me@275	147 description = ''
me@275	148 source = ''
me@275	149
me@275	150 def __init__(self, sequence=[], name=None, description=None, source=None):
me@275	151 super(Sequence, self).__init__(sequence)
me@275	152 if hasattr(sequence, 'name'):
me@275	153 vars(self).update(vars(sequence))
me@275	154 if name:
me@275	155 self.name = name
me@275	156 if description:
me@275	157 self.description = description
me@275	158 if source:
me@275	159 self.source = source
me@270	160
me@262	161 def __str__(self):
me@275	162 """Returns sequence in one-letter code."""
me@275	163 return ''.join(monomer.code1 for monomer in self)
me@270	164
me@273	165 @classmethod
me@273	166 def from_string(cls, string, name='', description=''):
me@273	167 """Create sequences from string of one-letter codes."""
me@273	168 monomer = cls.monomer_type.from_code1
me@273	169 monomers = [monomer(letter) for letter in string]
me@273	170 return cls(monomers, name, description)
me@262	171
me@275	172 @classmethod
me@275	173 def file_slice(cls, file, n_from, n_to, fasta_name='', name='', description=''):
me@270	174 """ Build and return sequence, consisting of part of sequence from file
me@262	175
me@262	176 Does not control gaps
me@262	177 """
me@262	178 inside = False
me@262	179 number_used = 0
me@262	180 s = ''
me@262	181 for line in file:
me@262	182 line = line.split()
me@262	183 if not inside:
me@262	184 if line.startswith('>%s' % fasta_name):
me@262	185 inside = True
me@262	186 else:
me@262	187 n = len(line)
me@262	188 s += line[(n_from - number_user):(n_to - number_user)]
me@275	189 return cls.from_string(s, name, description)
bnagaev@243	190
bnagaev@249	191 class Alignment(dict):
bnagaev@249	192 """ Alignment
me@270	193
bnagaev@249	194 {<Sequence object>:[<Monomer object>,None,<Monomer object>]}
me@266	195
bnagaev@249	196 keys are the Sequence objects, values are the lists, which
bnagaev@249	197 contain monomers of those sequences or None for gaps in the
bnagaev@249	198 corresponding sequence of alignment
bnagaev@249	199 """
bnagaev@249	200 # _sequences -- list of Sequence objects. Sequences don't contain gaps
bnagaev@249	201 # - see sequence.py module
bnagaev@249	202
bnagaev@249	203 def __init__(self, *args):
bnagaev@249	204 """overloaded constructor
bnagaev@249	205
bnagaev@249	206 Alignment() -> new empty Alignment
bnagaev@249	207 Alignment(sequences, body) -> new Alignment with sequences and
bnagaev@249	208 body initialized from arguments
bnagaev@249	209 Alignment(fasta_file) -> new Alignment, read body and sequences
me@270	210 from fasta file
bnagaev@249	211
bnagaev@249	212 """
bnagaev@249	213 if len(args)>1:#overloaded constructor
bnagaev@249	214 self.sequences=args[0]
bnagaev@249	215 self.body=args[1]
bnagaev@249	216 elif len(args)==0:
bnagaev@249	217 self.sequences=[]
bnagaev@249	218 self.body={}
bnagaev@249	219 else:
bnagaev@249	220 self.sequences, self.body = Alignment.from_fasta(args[0])
bnagaev@249	221
bnagaev@249	222 def length(self):
bnagaev@249	223 """ Returns width, ie length of each sequence with gaps """
bnagaev@249	224 return max([len(line) for line in self.body.values()])
bnagaev@249	225
bnagaev@249	226 def height(self):
bnagaev@249	227 """ The number of sequences in alignment (it's thickness). """
bnagaev@249	228 return len(self.body)
bnagaev@249	229
bnagaev@249	230 def identity(self):
bnagaev@249	231 """ Calculate the identity of alignment positions for colouring.
bnagaev@249	232
bnagaev@249	233 For every (row, column) in alignment the percentage of the exactly
bnagaev@249	234 same residue in the same column in the alignment is calculated.
me@270	235 The data structure is just like the Alignment.body, but istead of
bnagaev@249	236 monomers it contains float percentages.
bnagaev@249	237 """
bnagaev@249	238 # Oh, God, that's awful! Absolutely not understandable.
bnagaev@249	239 # First, calculate percentages of amino acids in every column
bnagaev@249	240 contribution = 1.0 / len(self.sequences)
bnagaev@249	241 all_columns = []
bnagaev@249	242 for position in range(len(self)):
bnagaev@249	243 column_percentage = {}
bnagaev@249	244 for seq in self.body:
bnagaev@249	245 if self.body[seq][position] is not None:
bnagaev@249	246 aa = self.body[seq][position].code
bnagaev@249	247 else:
bnagaev@249	248 aa = None
bnagaev@249	249 if aa in allpy.data.amino_acids:
bnagaev@249	250 if aa in column_percentage.keys():
bnagaev@249	251 column_percentage[aa] += contribution
bnagaev@249	252 else:
bnagaev@249	253 column_percentage[aa] = contribution
bnagaev@249	254 all_columns.append(column_percentage)
bnagaev@249	255 # Second, map these percentages onto the alignment
bnagaev@249	256 self.identity_percentages = {}
bnagaev@249	257 for seq in self.sequences:
bnagaev@249	258 self.identity_percentages[seq] = []
bnagaev@249	259 for seq in self.identity_percentages:
bnagaev@249	260 line = self.identity_percentages[seq]
bnagaev@249	261 for position in range(len(self)):
bnagaev@249	262 if self.body[seq][position] is not None:
bnagaev@249	263 aa = self.body[seq][position].code
bnagaev@249	264 else:
bnagaev@249	265 aa = None
bnagaev@249	266 line.append(all_columns[position].get(aa))
bnagaev@249	267 return self.identity_percentages
bnagaev@249	268
bnagaev@249	269 @staticmethod
bnagaev@249	270 def from_fasta(file, monomer_kind=AminoAcidType):
bnagaev@249	271 """ Import data from fasta file
me@270	272
bnagaev@249	273 monomer_kind is class, inherited from MonomerType
me@270	274
bnagaev@249	275 >>> import alignment
me@270	276 >>> sequences,body=alignment.Alignment.from_fasta(open("test.fasta"))
bnagaev@249	277 """
bnagaev@249	278 import re
bnagaev@249	279
bnagaev@249	280 sequences = []
bnagaev@249	281 body = {}
bnagaev@249	282
bnagaev@249	283 raw_sequences = file.read().split(">")
bnagaev@249	284 if len(raw_sequences) <= 1:
bnagaev@249	285 raise Exception("Wrong format of fasta-file %s" % file.name)
me@270	286
bnagaev@249	287 raw_sequences = raw_sequences[1:] #ignore everything before the first >
bnagaev@249	288 for raw in raw_sequences:
bnagaev@249	289 parsed_raw_sequence = raw.split("\n")
bnagaev@249	290 parsed_raw_sequence = [s.strip() for s in parsed_raw_sequence]
bnagaev@249	291 name_and_description = parsed_raw_sequence[0]
bnagaev@249	292 name_and_description = name_and_description.split(" ",1)
bnagaev@249	293 if len(name_and_description) == 2:
bnagaev@249	294 name, description = name_and_description
me@270	295 elif len(name_and_description) == 1:
bnagaev@249	296 #if there is description
bnagaev@249	297 name = name_and_description[0]
bnagaev@249	298 description = ''
bnagaev@249	299 else:
bnagaev@249	300 raise Exception("Wrong name of sequence %(name)$ fasta-file %(file)s" % \
bnagaev@249	301 {'name': name, 'file': file.name})
me@270	302
bnagaev@249	303 if len(parsed_raw_sequence) <= 1:
bnagaev@249	304 raise Exception("Wrong format of sequence %(name)$ fasta-file %(file)s" % \
bnagaev@249	305 {'name': name, 'file': file.name})
bnagaev@249	306 string = ""
bnagaev@249	307 for piece in parsed_raw_sequence[1:]:
bnagaev@249	308 piece_without_whitespace_chars = re.sub("\s", "", piece)
bnagaev@249	309 string += piece_without_whitespace_chars
bnagaev@249	310 monomers = [] #convert into Monomer objects
bnagaev@249	311 body_list = [] #create the respective list in body dict
bnagaev@249	312 for current_monomer in string:
bnagaev@249	313 if current_monomer not in ["-", ".", "~"]:
bnagaev@249	314 monomers.append(monomer_kind.from_code1(current_monomer).instance())
bnagaev@249	315 body_list.append(monomers[-1])
bnagaev@249	316 else:
bnagaev@249	317 body_list.append(None)
bnagaev@249	318 s = sequence.Sequence(monomers, name, description)
bnagaev@249	319 sequences.append(s)
bnagaev@249	320 body[s] = body_list
bnagaev@249	321 return sequences, body
me@270	322
bnagaev@249	323 @staticmethod
bnagaev@249	324 def from_sequences(*sequences):
bnagaev@249	325 """ Constructs new alignment from sequences
me@270	326
me@270	327 Add None's to right end to make equal lengthes of alignment sequences
bnagaev@249	328 """
bnagaev@249	329 alignment = Alignment()
bnagaev@249	330 alignment.sequences = sequences
bnagaev@249	331 max_length = max(len(sequence) for sequence in sequences)
bnagaev@249	332 for sequence in sequences:
bnagaev@249	333 gaps_count = max_length - len(sequence)
bnagaev@249	334 alignment.body[sequence] = sequence.monomers + [None] * gaps_count
bnagaev@249	335 return alignment
me@270	336
bnagaev@249	337 def save_fasta(self, out_file, long_line=70, gap='-'):
bnagaev@249	338 """ Saves alignment to given file
me@270	339
bnagaev@249	340 Splits long lines to substrings of length=long_line
me@270	341 To prevent this, set long_line=None
bnagaev@249	342 """
bnagaev@249	343 block.Block(self).save_fasta(out_file, long_line=long_line, gap=gap)
me@270	344
bnagaev@249	345 def muscle_align(self):
bnagaev@249	346 """ Simple align ths alignment using sequences (muscle)
me@270	347
bnagaev@249	348 uses old Monomers and Sequences objects
bnagaev@249	349 """
bnagaev@249	350 tmp_file = NamedTemporaryFile(delete=False)
bnagaev@249	351 self.save_fasta(tmp_file)
bnagaev@249	352 tmp_file.close()
bnagaev@249	353 os.system("muscle -in %(tmp)s -out %(tmp)s" % {'tmp': tmp_file.name})
bnagaev@249	354 sequences, body = Alignment.from_fasta(open(tmp_file.name))
bnagaev@249	355 for sequence in self.sequences:
bnagaev@249	356 try:
bnagaev@249	357 new_sequence = [i for i in sequences if sequence==i][0]
bnagaev@249	358 except:
bnagaev@249	359 raise Exception("Align: Cann't find sequence %s in muscle output" % \
bnagaev@249	360 sequence.name)
bnagaev@249	361 old_monomers = iter(sequence.monomers)
bnagaev@249	362 self.body[sequence] = []
bnagaev@249	363 for monomer in body[new_sequence]:
bnagaev@249	364 if not monomer:
bnagaev@249	365 self.body[sequence].append(monomer)
bnagaev@249	366 else:
bnagaev@249	367 old_monomer = old_monomers.next()
bnagaev@249	368 if monomer != old_monomer:
bnagaev@249	369 raise Exception("Align: alignment errors")
bnagaev@249	370 self.body[sequence].append(old_monomer)
bnagaev@249	371 os.unlink(tmp_file.name)
me@270	372
bnagaev@249	373 def column(self, sequence=None, sequences=None, original=None):
bnagaev@249	374 """ returns list of columns of alignment
me@270	375
bnagaev@249	376 sequence or sequences:
bnagaev@249	377 if sequence is given, then column is (original_monomer, monomer)
me@270	378 if sequences is given, then column is (original_monomer, {sequence: monomer})
bnagaev@249	379 if both of them are given, it is an error
bnagaev@249	380 original (Sequence type):
bnagaev@249	381 if given, this filters only columns represented by original sequence
bnagaev@249	382 """
bnagaev@249	383 if sequence and sequences:
bnagaev@249	384 raise Exception("Wrong usage. read help")
bnagaev@249	385 indexes = dict([(v, k) for( k, v) in enumerate(self.sequences)])
bnagaev@249	386 alignment = self.body.items()
bnagaev@249	387 alignment.sort(key=lambda i: indexes[i[0]])
bnagaev@249	388 alignment = [monomers for seq, monomers in alignment]
bnagaev@249	389 for column in zip(*alignment):
bnagaev@249	390 if not original or column[indexes[original]]:
bnagaev@249	391 if sequence:
bnagaev@249	392 yield (column[indexes[original]], column[indexes[sequence]])
bnagaev@249	393 else:
me@270	394 yield (column[indexes[original]],
bnagaev@249	395 dict([(s, column[indexes[s]]) for s in sequences]))
me@270	396
bnagaev@249	397 def secstr(self, sequence, pdb_chain, gap='-'):
bnagaev@249	398 """ Returns string representing secondary structure """
bnagaev@249	399 return ''.join([
me@270	400 (sequence.pdb_secstr[pdb_chain][m] if sequence.secstr_has(pdb_chain, m) else gap)
bnagaev@249	401 for m in self.body[sequence]])
bnagaev@249	402
bnagaev@249	403 class Block(object):
me@261	404 """ Block of alignment
me@270	405
me@261	406 Mandatory data:
me@266	407
me@261	408 * self.alignment -- alignment object, which the block belongs to
me@261	409 * self.sequences - set of sequence objects that contain monomers
me@261	410 and/or gaps, that constitute the block
me@261	411 * self.positions -- list of positions of the alignment.body that
me@261	412 are included in the block; position[i+1] is always to the right from position[i]
me@270	413
me@261	414 Don't change self.sequences -- it may be a link to other block.sequences
me@270	415
me@261	416 How to create a new block:
me@261	417 >>> import alignment
me@261	418 >>> import block
me@261	419 >>> proj = alignment.Alignment(open("test.fasta"))
me@261	420 >>> block1 = block.Block(proj)
me@261	421 """
me@270	422
me@261	423 def __init__(self, alignment, sequences=None, positions=None):
me@261	424 """ Builds new block from alignment
me@270	425
me@261	426 if sequences==None, all sequences are used
me@261	427 if positions==None, all positions are used
me@261	428 """
me@261	429 if sequences == None:
me@261	430 sequences = set(alignment.sequences) # copy
me@261	431 if positions == None:
me@261	432 positions = range(len(alignment))
me@261	433 self.alignment = alignment
me@261	434 self.sequences = sequences
me@261	435 self.positions = positions
me@270	436
me@261	437 def save_fasta(self, out_file, long_line=70, gap='-'):
me@270	438 """ Saves alignment to given file in fasta-format
me@270	439
me@261	440 No changes in the names, descriptions or order of the sequences
me@261	441 are made.
me@261	442 """
me@261	443 for sequence in self.sequences:
me@261	444 alignment_monomers = self.alignment.body[sequence]
me@261	445 block_monomers = [alignment_monomers[i] for i in self.positions]
me@261	446 string = ''.join([m.type.code1 if m else '-' for m in block_monomers])
me@261	447 save_fasta(out_file, string, sequence.name, sequence.description, long_line)
me@270	448
me@270	449 def geometrical_cores(self, max_delta=config.delta,
me@270	450 timeout=config.timeout, minsize=config.minsize,
me@261	451 ac_new_atoms=config.ac_new_atoms,
me@261	452 ac_count=config.ac_count):
me@261	453 """ Returns length-sorted list of blocks, representing GCs
me@270	454
me@261	455 max_delta -- threshold of distance spreading
me@261	456 timeout -- Bron-Kerbosh timeout (then fast O(n ln n) algorithm)
me@261	457 minsize -- min size of each core
me@261	458 ac_new_atoms -- min part or new atoms in new alternative core
me@261	459 current GC is compared with each of already selected GCs
me@261	460 if difference is less then ac_new_atoms, current GC is skipped
me@261	461 difference = part of new atoms in current core
me@261	462 ac_count -- max number of cores (including main core)
me@261	463 -1 means infinity
me@261	464 If more than one pdb chain for some sequence provided, consider all of them
me@270	465 cost is calculated as 1 / (delta + 1)
me@261	466 delta in [0, +inf) => cost in (0, 1]
me@261	467 """
me@261	468 nodes = self.positions
me@261	469 lines = {}
me@261	470 for i in self.positions:
me@261	471 for j in self.positions:
me@261	472 if i < j:
me@261	473 distances = []
me@261	474 for sequence in self.sequences:
me@261	475 for chain in sequence.pdb_chains:
me@261	476 m1 = self.alignment.body[sequence][i]
me@261	477 m2 = self.alignment.body[sequence][j]
me@261	478 if m1 and m2:
me@261	479 r1 = sequence.pdb_residues[chain][m1]
me@261	480 r2 = sequence.pdb_residues[chain][m2]
me@261	481 ca1 = r1['CA']
me@261	482 ca2 = r2['CA']
me@261	483 d = ca1 - ca2 # Bio.PDB feature
me@261	484 distances.append(d)
me@261	485 if len(distances) >= 2:
me@261	486 delta = max(distances) - min(distances)
me@261	487 if delta <= max_delta:
me@261	488 lines[Graph.line(i, j)] = 1.0 / (1.0 + max_delta)
me@261	489 graph = Graph(nodes, lines)
me@261	490 cliques = graph.cliques(timeout=timeout, minsize=minsize)
me@261	491 GCs = []
me@261	492 for clique in cliques:
me@261	493 for GC in GCs:
me@261	494 if len(clique - set(GC.positions)) < ac_new_atoms * len(clique):
me@261	495 break
me@261	496 else:
me@261	497 GCs.append(Block(self.alignment, self.sequences, clique))
me@261	498 if ac_count != -1 and len(GCs) >= ac_count:
me@261	499 break
me@261	500 return GCs
me@270	501
me@261	502 def xstring(self, x='X', gap='-'):
me@261	503 """ Returns string consisting of gap chars and chars x at self.positions
me@270	504
me@261	505 Length of returning string = length of alignment
me@261	506 """
me@261	507 monomers = [False] * len(self.alignment)
me@261	508 for i in self.positions:
me@261	509 monomers[i] = True
me@261	510 return ''.join([x if m else gap for m in monomers])
me@270	511
me@261	512 def save_xstring(self, out_file, name, description='', x='X', gap='-', long_line=70):
me@261	513 """ Save xstring and name in fasta format """
me@261	514 save_fasta(out_file, self.xstring(x=x, gap=gap), name, description, long_line)
me@270	515
me@261	516 def monomers(self, sequence):
me@261	517 """ Iterates monomers of this sequence from this block """
me@261	518 alignment_sequence = self.alignment.body[sequence]
me@261	519 return (alignment_sequence[i] for i in self.positions)
me@270	520
me@261	521 def ca_atoms(self, sequence, pdb_chain):
me@261	522 """ Iterates Ca-atom of monomers of this sequence from this block """
me@261	523 return (sequence.pdb_residues[pdb_chain][monomer] for monomer in self.monomers())
me@270	524
me@261	525 def sequences_chains(self):
me@261	526 """ Iterates pairs (sequence, chain) """
me@261	527 for sequence in self.alignment.sequences:
me@261	528 if sequence in self.sequences:
me@261	529 for chain in sequence.pdb_chains:
me@261	530 yield (sequence, chain)
me@270	531
me@261	532 def superimpose(self):
me@261	533 """ Superimpose all pdb_chains in this block """
me@261	534 sequences_chains = list(self.sequences_chains())
me@261	535 if len(sequences_chains) >= 1:
me@261	536 sup = Superimposer()
me@261	537 fixed_sequence, fixed_chain = sequences_chains.pop()
me@261	538 fixed_atoms = self.ca_atoms(fixed_sequence, fixed_chain)
me@261	539 for sequence, chain in sequences_chains:
me@261	540 moving_atoms = self.ca_atoms(sequence, chain)
me@261	541 sup.set_atoms(fixed_atoms, moving_atoms)
me@261	542 # Apply rotation/translation to the moving atoms
me@261	543 sup.apply(moving_atoms)
me@270	544
me@261	545 def pdb_save(self, out_file):
me@270	546 """ Save all sequences
me@270	547
me@261	548 Returns {(sequence, chain): CHAIN}
me@261	549 CHAIN is chain letter in new file
me@261	550 """
me@261	551 tmp_file = NamedTemporaryFile(delete=False)
me@261	552 tmp_file.close()
me@270	553
me@261	554 for sequence, chain in self.sequences_chains():
me@261	555 sequence.pdb_save(tmp_file.name, chain)
me@261	556 # TODO: read from tmp_file.name
me@261	557 # change CHAIN
me@261	558 # add to out_file
me@270	559
me@261	560 os.unlink(NamedTemporaryFile)
bnagaev@239	561
me@260	562 # vim: set ts=4 sts=4 sw=4 et: