Package Bio :: Package PopGen :: Package SimCoal :: Module Template
[hide private]
[frames] | no frames]

Source Code for Module Bio.PopGen.SimCoal.Template

  1  # Copyright 2007 by Tiago Antao.  All rights reserved. 
  2  # This code is part of the Biopython distribution and governed by its 
  3  # license.  Please see the LICENSE file that should have been included 
  4  # as part of this package. 
  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   
14 -def exec_template(template):
15 executed_template = template 16 match = re.search('!!!(.*?)!!!', executed_template, re.MULTILINE) 17 #while len(match.groups())>0: 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 #match = patt.matcher(String(executed_template)) 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 #f_name += '_' + str(total_size) 33 for tup in curr_values: 34 name, val = tup 35 f_name += '_' + str(val) 36 #reg = re.compile('\?' + name, re.MULTILINE) 37 #template = re.sub(reg, str(val), template) 38 template = template.replace('?'+name, str(val)) 39 f = open(f_name + '.par', 'w') 40 #executed_template = template 41 executed_template = exec_template(template) 42 clean_template = executed_template.replace('\r\n','\n').replace('\n\n','\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
63 -def get_xy_from_matrix(x_max, y_max, pos):
64 y = (pos-1) / x_max 65 x = (pos-1) % x_max 66 return x, y
67
68 -def get_step_2d(x_max, y_max, x, y, mig):
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 -def generate_ssm2d_mat(x_max, y_max, mig):
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
85 -def generate_island_mat(total_size, mig):
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
96 -def generate_null_mat(total_size):
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
104 -def generate_join_events(t, total_size, join_size, orig_size):
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
111 -def no_processor(in_string):
112 return in_string
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 #def prepare_dir(): 120 # try: 121 # mkdir(sep.join([Config.dataDir, 'SimCoal'])) #Should exist, but... 122 # except OSError: 123 # pass #Its ok if already exists 124 # try: 125 # mkdir(sep.join([Config.dataDir, 'SimCoal', 'runs'])) 126 # except OSError: 127 # pass #Its ok if already exists 128 129 130 #sep is because of jython
131 -def generate_model(par_stream, out_prefix, params, 132 specific_processor = no_processor, out_dir = '.'):
133 #prepare_dir() 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
139 -def get_demography_template(stream, model, tp_dir = None):
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 #Internal Template 152 f = open(sep.join([builtin_tpl_dir, model + '.par']), 'r') 153 else: 154 #External template 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
162 -def _gen_loci(stream, loci):
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
171 -def get_chr_template(stream, chrs):
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
200 -def generate_simcoal_from_template(model, chrs, params, out_dir = '.', tp_dir=None):
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 #par_stream = open(out_dir + sep + 'tmp.par', 'r') 215 #print par_stream.read() 216 #par_stream.close() 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