1
2
3
4
5
6 from math import sqrt
7 from sys import argv,exit
8 from os import sep, mkdir
9 import re
10
11 from Bio.PopGen.SimCoal import builtin_tpl_dir
12
13
15 executed_template = template
16 match = re.search('!!!(.*?)!!!', executed_template, re.MULTILINE)
17
18 while match :
19 exec_result = str(eval(match.groups()[0]))
20 executed_template = executed_template.replace(
21 '!!!' + match.groups()[0] + '!!!',
22 exec_result, 1)
23 match = re.search('!!!(.*?)!!!', executed_template, re.MULTILINE)
24
25 return executed_template
26
27
28 -def process_para(in_string, out_file_prefix, para_list, curr_values):
29 if (para_list == []):
30 template = in_string
31 f_name = out_file_prefix
32
33 for tup in curr_values:
34 name, val = tup
35 f_name += '_' + str(val)
36
37
38 template = template.replace('?'+name, str(val))
39 f = open(f_name + '.par', 'w')
40
41 executed_template = exec_template(template)
42 clean_template = executed_template.replace('\r\n','\n').replace('\n\n','\n').replace('\n','\r\n')
43 f.write(clean_template)
44 f.close()
45 return [f_name]
46 else:
47 name, rng = para_list[0]
48 fnames = []
49 for val in rng:
50 new_values = [(name, val)]
51 new_values.extend(curr_values)
52 more_names = process_para(in_string, out_file_prefix, para_list[1:], new_values)
53 fnames.extend(more_names)
54 return fnames
55
56
57 -def dupe(motif, times):
58 ret_str = ''
59 for i in range(1, times + 1):
60 ret_str += motif + '\r\n'
61 return ret_str
62
64 y = (pos-1) / x_max
65 x = (pos-1) % x_max
66 return x, y
67
69 my_x, my_y = get_xy_from_matrix(x_max, y_max, y)
70 other_x, other_y = get_xy_from_matrix(x_max, y_max, x)
71
72 if (my_x-other_x)**2 + (my_y-other_y)**2 == 1:
73 return str(mig) + ' '
74 else:
75 return '0 '
76
78 mig_mat = ''
79 for x in range(1, x_max*y_max + 1):
80 for y in range(1, x_max*y_max + 1):
81 mig_mat += get_step_2d(x_max, y_max, x, y, mig)
82 mig_mat += "\r\n"
83 return mig_mat
84
86 mig_mat = ''
87 for x in range(1, total_size + 1):
88 for y in range(1, total_size + 1):
89 if (x == y):
90 mig_mat += '0 '
91 else:
92 mig_mat += '!!!' + str(mig) + '!!! '
93 mig_mat += "\r\n"
94 return mig_mat
95
97 null_mat = ''
98 for x in range(1, total_size + 1):
99 for y in range(1, total_size + 1):
100 null_mat += '0 '
101 null_mat += '\r\n'
102 return null_mat
103
105 events = ''
106 for i in range(1, total_size-1):
107 events += str(t) + ' ' + str(i) + ' 0 1 1 0 1\r\n'
108 events += str(t) + ' ' + str(total_size-1) + ' 0 1 ' + str(1.0*total_size*join_size/orig_size) + ' 0 1\r\n'
109 return events
110
113
114 -def process_text(in_string, out_file_prefix, para_list, curr_values,
115 specific_processor):
116 text = specific_processor(in_string)
117 return process_para(text, out_file_prefix, para_list, [])
118
119
120
121
122
123
124
125
126
127
128
129
130
133
134 text = par_stream.read()
135 out_file_prefix = sep.join([out_dir, out_prefix])
136 return process_text(text, out_file_prefix, params, [], specific_processor)
137
138
140 '''
141 Gets a demograpy template.
142
143 Most probably this model needs to be sent to GenCases.
144
145 stream - Writable stream.
146 param - Template file.
147 tp_dir - Directory where to find the template, if None
148 use an internal template
149 '''
150 if tp_dir == None:
151
152 f = open(sep.join([builtin_tpl_dir, model + '.par']), 'r')
153 else:
154
155 f = open(sep.join([tp_dir, model + '.par']), 'r')
156 l = f.readline()
157 while l<>'':
158 stream.write(l)
159 l = f.readline()
160 f.close()
161
163 stream.write('//Number of contiguous linkage blocks in chromosome\n')
164 stream.write(str(len(loci)) + '\n')
165 stream.write('//Per Block: Data type, No. of loci, Recombination rate to the right-side locus, plus optional parameters\n')
166 for locus in loci:
167 stream.write(' '.join([locus[0]] +
168 map(lambda x: str(x), list(locus[1])
169 )) + '\n')
170
172 '''
173 Writes a Simcoal2 loci template part.
174
175 stream - Writable stream.
176 chr - Chromosome list.
177
178 Current loci list:
179 [(chr_repeats,[(marker, (params))])]
180 chr_repeats --> Number of chromosome repeats
181 marker --> 'SNP', 'DNA', 'RFLP', 'MICROSAT'
182 params --> Simcoal2 parameters for markers (list of floats
183 or ints - if to be processed by generate_model)
184 '''
185 num_chrs = reduce(lambda x, y: x + y[0], chrs, 0)
186 stream.write('//Number of independent (unlinked) chromosomes, and "chromosome structure" flag: 0 for identical structure across chromosomes, and 1 for different structures on different chromosomes.\n')
187 if len(chrs) > 1 or num_chrs == 1:
188 stream.write(str(num_chrs) + ' 1\n')
189 else:
190 stream.write(str(num_chrs) + ' 0\n')
191 for chr in chrs:
192 repeats = chr[0]
193 loci = chr[1]
194 if len(chrs) == 1:
195 _gen_loci(stream, loci)
196 else:
197 for i in range(repeats):
198 _gen_loci(stream, loci)
199
201 '''
202 Writes a complete SimCoal2 template file.
203
204 This joins together get_demography_template and get_chr_template,
205 which are feed into generate_model
206 Please check the three functions for parameters (model from
207 get_demography_template, chrs from get_chr_template and
208 params from generate_model).
209 '''
210 stream = open(out_dir + sep + 'tmp.par', 'w')
211 get_demography_template(stream, model, tp_dir)
212 get_chr_template(stream, chrs)
213 stream.close()
214
215
216
217 par_stream = open(out_dir + sep + 'tmp.par', 'r')
218 generate_model(par_stream, model, params, out_dir = out_dir)
219 par_stream.close()
220