Back to index

python-biopython  1.60
Template.py
Go to the documentation of this file.
00001 # Copyright 2007 by Tiago Antao.  All rights reserved.
00002 # This code is part of the Biopython distribution and governed by its
00003 # license.  Please see the LICENSE file that should have been included
00004 # as part of this package.
00005 
00006 from math import sqrt
00007 from sys import argv,exit
00008 from os import sep, mkdir
00009 import re
00010 
00011 from Bio.PopGen.SimCoal import builtin_tpl_dir
00012 
00013 
00014 def exec_template(template):
00015     executed_template = template
00016     match = re.search('!!!(.*?)!!!', executed_template, re.MULTILINE)
00017     #while len(match.groups())>0:
00018     while match:
00019         exec_result = str(eval(match.groups()[0]))
00020         executed_template = executed_template.replace(
00021                '!!!' + match.groups()[0] + '!!!',
00022                exec_result, 1)
00023         match = re.search('!!!(.*?)!!!', executed_template, re.MULTILINE)
00024         #match = patt.matcher(String(executed_template))
00025     return executed_template
00026     
00027 
00028 def process_para(in_string, out_file_prefix, para_list, curr_values):
00029     if (para_list == []):
00030         template = in_string
00031         f_name = out_file_prefix
00032         #f_name += '_' + str(total_size)
00033         for tup in curr_values:
00034             name, val = tup
00035             f_name += '_' + str(val) 
00036             #reg = re.compile('\?' + name, re.MULTILINE)
00037             #template = re.sub(reg, str(val), template)
00038             template = template.replace('?'+name, str(val))
00039         f = open(f_name + '.par', 'w')
00040         #executed_template = template
00041         executed_template = exec_template(template)
00042         clean_template =  executed_template.replace('\r\n','\n').replace('\n\n','\n')
00043         f.write(clean_template)
00044         f.close()
00045         return [f_name]
00046     else:
00047         name, rng = para_list[0]
00048         fnames = []
00049         for val in rng:
00050             new_values = [(name, val)]
00051             new_values.extend(curr_values)
00052             more_names = process_para(in_string, out_file_prefix, para_list[1:], new_values)
00053             fnames.extend(more_names)
00054         return fnames
00055 
00056 
00057 def dupe(motif, times):
00058     ret_str = ''
00059     for i in range(1, times + 1):
00060         ret_str += motif + '\r\n'
00061     return ret_str
00062 
00063 def get_xy_from_matrix(x_max, y_max, pos):
00064     y = (pos-1) / x_max
00065     x = (pos-1) % x_max
00066     return x, y
00067 
00068 def get_step_2d(x_max, y_max, x, y, mig):
00069     my_x,    my_y    = get_xy_from_matrix(x_max, y_max, y)
00070     other_x, other_y = get_xy_from_matrix(x_max, y_max, x)
00071         
00072     if (my_x-other_x)**2 + (my_y-other_y)**2 == 1:
00073         return str(mig) + ' '
00074     else:
00075         return '0 '
00076 
00077 def generate_ssm2d_mat(x_max, y_max, mig):
00078     mig_mat = ''
00079     for x in range(1, x_max*y_max + 1):
00080         for y in range(1, x_max*y_max + 1):
00081             mig_mat += get_step_2d(x_max, y_max, x, y, mig)
00082         mig_mat += "\r\n"
00083     return mig_mat
00084 
00085 def generate_island_mat(total_size, mig):
00086     mig_mat = ''
00087     for x in range(1, total_size + 1):
00088         for y in range(1, total_size + 1):
00089             if (x == y):
00090                 mig_mat += '0 '
00091             else:
00092                 mig_mat += '!!!' + str(mig) + '!!! '
00093         mig_mat += "\r\n"
00094     return mig_mat
00095 
00096 def generate_null_mat(total_size):
00097     null_mat = ''
00098     for x in range(1, total_size + 1):
00099         for y in range(1, total_size + 1):
00100             null_mat += '0 '
00101         null_mat += '\r\n'
00102     return null_mat
00103 
00104 def generate_join_events(t, total_size, join_size, orig_size):
00105     events = ''
00106     for i in range(1, total_size-1):
00107         events += str(t) + ' ' + str(i) + ' 0 1 1 0 1\r\n'
00108     events += str(t) + ' ' + str(total_size-1) + ' 0 1 ' + str(1.0*total_size*join_size/orig_size) + ' 0 1\r\n'
00109     return events
00110 
00111 def no_processor(in_string):
00112     return in_string
00113 
00114 def process_text(in_string, out_file_prefix, para_list, curr_values,
00115                  specific_processor):
00116     text = specific_processor(in_string)
00117     return process_para(text, out_file_prefix, para_list, [])
00118 
00119 #def prepare_dir():
00120 #    try:
00121 #        mkdir(sep.join([Config.dataDir, 'SimCoal'])) #Should exist, but...
00122 #    except OSError:
00123 #        pass #Its ok if already exists
00124 #    try:
00125 #        mkdir(sep.join([Config.dataDir, 'SimCoal', 'runs']))
00126 #    except OSError:
00127 #        pass #Its ok if already exists
00128     
00129 
00130 #sep is because of jython
00131 def generate_model(par_stream, out_prefix, params,
00132     specific_processor = no_processor, out_dir = '.'):
00133     #prepare_dir()
00134     text = par_stream.read()
00135     out_file_prefix = sep.join([out_dir, out_prefix])
00136     return process_text(text, out_file_prefix, params, [], specific_processor)
00137 
00138 
00139 def get_demography_template(stream, model, tp_dir = None):
00140     '''
00141         Gets a demograpy template.
00142  
00143         Most probably this model needs to be sent to GenCases.
00144  
00145         stream - Writable stream.
00146         param  - Template file.
00147         tp_dir - Directory where to find the template, if None
00148                  use an internal template
00149     '''
00150     if tp_dir == None:
00151         #Internal Template
00152         f = open(sep.join([builtin_tpl_dir, model + '.par']), 'r')
00153     else:
00154         #External template
00155         f = open(sep.join([tp_dir, model + '.par']), 'r')
00156     l = f.readline()
00157     while l!='':
00158         stream.write(l)
00159         l = f.readline()
00160     f.close()
00161 
00162 def _gen_loci(stream, loci):
00163     stream.write('//Number of contiguous linkage blocks in chromosome\n')
00164     stream.write(str(len(loci)) + '\n')
00165     stream.write('//Per Block: Data type, No. of loci, Recombination rate to the right-side locus, plus optional parameters\n')
00166     for locus in loci:
00167         stream.write(' '.join([locus[0]] +
00168             map(lambda x: str(x), list(locus[1])
00169         )) + '\n')
00170 
00171 def get_chr_template(stream, chrs):
00172     '''
00173         Writes a Simcoal2 loci template part.
00174 
00175         stream - Writable stream.
00176         chr    - Chromosome list.
00177 
00178         Current loci list:
00179           [(chr_repeats,[(marker, (params))])]
00180           chr_repeats --> Number of chromosome repeats
00181           marker  --> 'SNP', 'DNA', 'RFLP', 'MICROSAT'
00182           params  --> Simcoal2 parameters for markers (list of floats
00183             or ints - if to be processed by generate_model)
00184     '''
00185     num_chrs = reduce(lambda x, y: x + y[0], chrs, 0)
00186     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')
00187     if len(chrs) > 1 or num_chrs == 1:
00188         stream.write(str(num_chrs) + ' 1\n')
00189     else:
00190         stream.write(str(num_chrs) + ' 0\n')
00191     for chr in chrs:
00192         repeats = chr[0]
00193         loci = chr[1]
00194         if len(chrs) == 1:
00195             _gen_loci(stream, loci)
00196         else:
00197             for i in range(repeats):
00198                 _gen_loci(stream, loci)
00199 
00200 def generate_simcoal_from_template(model, chrs, params, out_dir = '.', tp_dir=None):
00201     '''
00202        Writes a complete SimCoal2 template file.
00203 
00204        This joins together get_demography_template and get_chr_template,
00205        which are feed into generate_model
00206        Please check the three functions for parameters (model from
00207          get_demography_template, chrs from get_chr_template and
00208          params from generate_model).
00209     '''
00210     stream = open(out_dir + sep + 'tmp.par', 'w')
00211     get_demography_template(stream, model, tp_dir)
00212     get_chr_template(stream, chrs)
00213     stream.close()
00214     #par_stream = open(out_dir + sep + 'tmp.par', 'r')
00215     #print par_stream.read()
00216     #par_stream.close()
00217     par_stream = open(out_dir + sep + 'tmp.par', 'r')
00218     generate_model(par_stream, model, params, out_dir = out_dir)
00219     par_stream.close()
00220 
00221