Документ взят из кэша поисковой машины. Адрес оригинального документа : http://kodomo.fbb.msu.ru/hg/allpy/file/db9d116e979f/lib/block.py
Дата изменения: Unknown
Дата индексирования: Sun Feb 3 17:56:42 2013
Кодировка:
allpy: db9d116e979f lib/block.py

allpy

view lib/block.py @ 147:db9d116e979f

find changesets by author, revision, files, or words in the commit message
documentation improvements
author boris <bnagaev@gmail.com>
date Sun, 24 Oct 2010 23:12:15 +0400
parents 8f1d8ece31af
children bc4c291cc2b0
line source
1 #!usr/bin/python
2
3 import sys
4
5 import project
6 import sequence
7 import monomer
8 import config
9 from graph import Graph
10 from Bio.PDB import Superimposer
11
12 class Block(object):
13 """ Block of alignment
14
15 Mandatory data:
16 * self.project -- project object, which the block belongs to
17 * self.sequences - set of sequence objects that contain monomers
18 and/or gaps, that constitute the block
19 * self.positions -- sorted list of positions of the project.alignment that
20 are included in the block
21
22 Don't change self.sequences -- it may be a link to other block.sequences
23
24 How to create a new block:
25 >>> import project
26 >>> import block
27 >>> proj = project.Project(open("test.fasta"))
28 >>> block1 = block.Block(proj)
29 """
30
31 def __init__(self, project, sequences=None, positions=None):
32 """ Builds new block from project
33
34 if sequences==None, all sequences are used
35 if positions==None, all positions are used
36 """
37 if sequences == None:
38 sequences = set(project.sequences) # copy
39 if positions == None:
40 positions = range(len(project))
41 self.project = project
42 self.sequences = sequences
43 self.positions = positions
44
45 def save_fasta(self, out_file, long_line=60, gap='-'):
46 """
47 Saves alignment to given file in fasta-format
48 Splits long lines to substrings of length=long_line
49 To prevent this, set long_line=None
50
51 No changes in the names, descriptions or order of the sequences
52 are made.
53 """
54 for sequence in self.sequences:
55 out_file.write(">%(name)s %(description)s \n" % sequence.__dict__)
56 alignment_monomers = self.project.alignment[sequence]
57 block_monomers = [alignment_monomers[i] for i in self.positions]
58 string = ''.join([m.type.code1 if m else '-' for m in block_monomers])
59 if long_line:
60 for i in range(0, len(string) // long_line + 1):
61 out_file.write("%s \n" % string[i*long_line : i*long_line + long_line])
62 else:
63 out_file.write("%s \n" % string)
64
65 def geometrical_cores(self, max_delta=config.delta,
66 timeout=config.timeout, minsize=config.minsize,
67 ac_new_atoms=config.ac_new_atoms,
68 ac_count=config.ac_count):
69 """
70 returns length-sorted list of blocks, representing GCs
71
72 max_delta -- threshold of distance spreading
73 timeout -- Bron-Kerbosh timeout (then fast O(n ln n) algorithm)
74 minsize -- min size of each core
75 ac_new_atoms -- min part or new atoms in new alternative core
76 current GC is compared with each of already selected GCs
77 if difference is less then ac_new_atoms, current GC is skipped
78 difference = part of new atoms in current core
79 ac_count -- max number of cores (including main core)
80 -1 means infinity
81 If more than one pdb chain for some sequence provided, consider all of them
82 cost is calculated as 1 / (delta + 1)
83 delta in [0, +inf) => cost in (0, 1]
84 """
85 nodes = self.positions
86 lines = {}
87 for i in self.positions:
88 for j in self.positions:
89 if i < j:
90 distances = []
91 for sequence in self.sequences:
92 for chain in sequence.pdb_chains:
93 m1 = self.project.alignment[sequence][i]
94 m2 = self.project.alignment[sequence][j]
95 if m1 and m2:
96 ca1 = m1.pdb_residues[chain]['CA']
97 ca2 = m2.pdb_residues[chain]['CA']
98 d = ca1 - ca2 # Bio.PDB feature
99 distances.append(d)
100 if len(distances) >= 2:
101 delta = max(distances) - min(distances)
102 if delta <= max_delta:
103 lines[Graph.line(i, j)] = 1.0 / (1.0 + max_delta)
104 graph = Graph(nodes, lines)
105 cliques = graph.cliques(timeout=timeout, minsize=minsize)
106 GCs = []
107 for clique in cliques:
108 for GC in GCs:
109 if len(clique - set(GC.positions)) < ac_new_atoms * len(clique):
110 break
111 else:
112 GCs.append(Block(self.project, self.sequences, clique))
113 if ac_count != -1 and len(GCs) >= ac_count:
114 break
115 return GCs
116
117 def xstring(self, x='X', gap='-'):
118 """
119 Returns string consisting of gap chars and chars x at self.positions
120 Length of returning string = length of project
121 """
122 monomers = [False] * len(self.project)
123 for i in self.positions:
124 monomers[i] = True
125 return ''.join([x if m else gap for m in monomers])
126
127 def save_xstring(self, out_file, name, description='', x='X', gap='-'):
128 """
129 Save xstring and name in fasta format
130 """
131 out_file.write(">%(name)s %(description)s \n" % \
132 {'name':name, 'description':description})