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