allpy: lib/block.py annotate

allpy

annotate lib/block.py @ 123:b7c3691c47bf

lib::block::xstring (rename it, please)

author	boris <bnagaev@gmail.com>
date	Sat, 23 Oct 2010 23:29:18 +0400
parents	48dc5974eead
children	9f96bc38bc3d

rev	line source
BurkovBA@0	1 #!usr/bin/python
BurkovBA@0	2
BurkovBA@0	3 import sys
BurkovBA@0	4
BurkovBA@0	5 import project
BurkovBA@0	6 import sequence
BurkovBA@0	7 import monomer
bnagaev@116	8 import config
bnagaev@121	9 from graph import Graph
BurkovBA@0	10
BurkovBA@0	11 class Block(object):
BurkovBA@0	12 """
BurkovBA@0	13 Mandatory data:
BurkovBA@0	14 * self.project -- project object, which the block belongs to
BurkovBA@1	15 * self.sequences - set of sequence objects that contain monomers
BurkovBA@0	16 and/or gaps, that constitute the block
bnagaev@115	17 * self.positions -- sorted list of positions of the project.alignment that
BurkovBA@1	18 are included in the block
bnagaev@116	19
BurkovBA@0	20 How to create a new block:
BurkovBA@0	21 >>> import project
BurkovBA@0	22 >>> import block
BurkovBA@0	23 >>> proj = project.Project(open("test.fasta"))
bnagaev@114	24 >>> block1 = block.Block(proj)
BurkovBA@0	25 """
BurkovBA@0	26
bnagaev@112	27 def __init__(self, project, sequences=None, positions=None):
bnagaev@112	28 """
bnagaev@112	29 Builds new block from project
bnagaev@112	30 if sequences==None, all sequences are used
bnagaev@112	31 if positions==None, all positions are used
bnagaev@112	32 """
bnagaev@112	33 if sequences == None:
bnagaev@112	34 sequences = project.sequences
bnagaev@112	35 if positions == None:
bnagaev@112	36 positions = range(len(project))
BurkovBA@73	37 self.project = project
BurkovBA@73	38 self.sequences = sequences
BurkovBA@73	39 self.positions = positions
BurkovBA@0	40
bnagaev@112	41 def save_fasta(self, out_file, long_line=60):
bnagaev@112	42 """
bnagaev@112	43 Saves alignment to given file in fasta-format
bnagaev@112	44 Splits long lines to substrings of length=long_line
bnagaev@112	45 To prevent this, set long_line=None
bnagaev@112	46
BurkovBA@0	47 No changes in the names, descriptions or order of the sequences
BurkovBA@0	48 are made.
BurkovBA@0	49 """
BurkovBA@0	50 for sequence in self.sequences:
bnagaev@112	51 out_file.write(">%(name)s %(description)s \n" % sequence.__dict__)
bnagaev@113	52 alignment_monomers = self.project.alignment[sequence]
bnagaev@115	53 block_monomers = [alignment_monomers[i] for i in self.positions]
bnagaev@113	54 string = ''.join([m.type.code1 if m else '-' for m in block_monomers])
bnagaev@112	55 if long_line:
bnagaev@112	56 for i in range(0, len(string) // long_line + 1):
bnagaev@112	57 out_file.write("%s \n" % string[ilong_line : ilong_line + long_line])
bnagaev@112	58 else:
bnagaev@112	59 out_file.write("%s \n" % string)
bnagaev@116	60
bnagaev@120	61 def geometrical_core(self, max_delta=config.delta,
bnagaev@120	62 timeout=config.timeout, minsize=config.minsize):
bnagaev@116	63 """
bnagaev@116	64 returns sorted list of positions, representing geometrical core
bnagaev@116	65 delta -- threshold of distance spreading
bnagaev@120	66
bnagaev@120	67 If more than one pdb chain for some sequence provided, consider all of them
bnagaev@120	68 cost is calculated as 1 / (delta + 1)
bnagaev@120	69 delta in [0, +inf) => cost in (0, 1]
bnagaev@116	70 """
bnagaev@117	71 nodes = self.positions
bnagaev@117	72 lines = {}
bnagaev@116	73 for i in self.positions:
bnagaev@116	74 for j in self.positions:
bnagaev@117	75 if i < j:
bnagaev@120	76 distances = []
bnagaev@120	77 for sequence in self.sequences:
bnagaev@120	78 for chain in sequence.pdb_chains:
bnagaev@120	79 m1 = self.project.alignment[sequence][i]
bnagaev@120	80 m2 = self.project.alignment[sequence][j]
bnagaev@122	81 if m1 and m2:
bnagaev@122	82 ca1 = m1.pdb_residues[chain]['CA']
bnagaev@122	83 ca2 = m2.pdb_residues[chain]['CA']
bnagaev@122	84 d = ca1 - ca2 # Bio.PDB feature
bnagaev@122	85 distances.append(d)
bnagaev@122	86 if len(distances) >= 2:
bnagaev@122	87 delta = max(distances) - min(distances)
bnagaev@122	88 if delta <= max_delta:
bnagaev@122	89 lines[Graph.line(i, j)] = 1.0 / (1.0 + max_delta)
bnagaev@120	90 graph = Graph(nodes, lines)
bnagaev@120	91 return graph.cliques(timeout=timeout, minsize=minsize)
bnagaev@123	92
bnagaev@123	93 def xstring(self, x):
bnagaev@123	94 """
bnagaev@123	95 Returns string consisting of '-' and chars x at self.positions
bnagaev@123	96 Length of returning string = length of project
bnagaev@123	97 """
bnagaev@123	98 monomers = [False] * len(self.project)
bnagaev@123	99 for i in self.positions:
bnagaev@123	100 monomers[i] = True
bnagaev@123	101 return ''.join([x if m else '-' for m in monomers])