allpy: db9d116e979f lib/block.py

allpy

view lib/block.py @ 147:db9d116e979f

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

3 import sys

5 import project

6 import sequence

7 import monomer

8 import config

9 from graph import Graph

10 from Bio.PDB import Superimposer

12 class Block(object):

13 """ Block of alignment

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

22 Don't change self.sequences -- it may be a link to other block.sequences

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 """

31 def __init__(self, project, sequences=None, positions=None):

32 """ Builds new block from project

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

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

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)

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

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})

133

134 out_file.write("%(xstring)s \n" % {'xstring':self.xstring(x=x, gap=gap)})

135

136 def monomers(self, sequence):

137 """ Iterates monomers of this sequence from this block """

138 alignment_sequence = self.project.alignment[sequence]

139 return (alignment_sequence[i] for i in self.positions)

140

141 def ca_atoms(self, sequence, pdb_chain):

142 """ Iterates Ca-atom of monomers of this sequence from this block """

143 return (monomer.pdb_residues[pdb_chain] for monomer in self.monomers())

144

145 def sequences_chains(self):

146 """ Iterates pairs (sequence, chain) """

147 for sequence in self.sequences:

148 for chain in sequence.pdb_chains:

149 yield (sequence, chain)

150

151 def superimpose(self):

152 """ Superimpose all pdb_chains in this block """

153 sequences_chains = []

154 if len(self.sequences) >= 1:

155 sup = Superimposer()

156 fixed_sequence, fixed_chain = sequences_chains.pop()

157 fixed_atoms = self.ca_atoms(fixed_sequence, fixed_chain)

158 for sequence, chain in sequences_chains:

159 moving_atoms = self.ca_atoms(sequence, chain)

160 sup.set_atoms(fixed_atoms, moving_atoms)

161 # Apply rotation/translation to the moving atoms

162 sup.apply(moving_atoms)

163

164 def pdb_save(self, out_file):

165 """ Superimpose all pdb_chains in this block """

166 tmp_file, tmp_filename = mkstemp()

167 os.close(tmp_file) # this is file descriptor, not normal file object.

168 sequence.pdb_save(tmp_filename, chain)

169

170

171