3 import os, string, sys, getopt, Args, optparse
5 record_delimiter = "\n---\n"
19 ###########################################
22 ###########################################
24 ###########################################
25 def set_values(self, win_size, step):
28 self.win_size = win_size
31 ###########################################
33 ###########################################
34 def open_streams(self, input_file, output_file):
36 self.in_stream = sys.stdin
39 self.in_stream = open(input_file, 'r')
44 self.out_stream = sys.stdout
47 self.out_stream = open(output_file, 'w')
50 ###########################################
52 ###########################################
53 def close_streams(self):
55 self.in_stream.close()
57 self.out_stream.close()
58 ###########################################
60 ###########################################
64 while not self.eo_buffer:
65 if self.in_stream.isatty():
73 tmp = self.in_stream.read(1000)
77 self.buffer = self.buffer + tmp
78 delim_index = self.buffer.find(record_delimiter)
80 rec = self.buffer[:delim_index]
81 self.buffer = self.buffer[delim_index + len(record_delimiter):]
84 ###########################################
86 ###########################################
87 def process_record(self, rec):
88 lines = rec.split("\n")
90 self.rec_dic[self.rec_num] = {}
98 self.rec_dic[self.rec_num][toks[0]] = toks[1]
100 score = string.atoi(toks[1])
102 chr = self.rec_dic[self.rec_num]["CHR"]
103 chr_beg = string.atoi(self.rec_dic[self.rec_num]["CHR_BEG"])
105 if not self.density_dic.has_key(chr):
106 self.density_dic[chr] = {}
108 k_cur = (chr_beg / step) * step
109 if chr_beg in range(k_cur, k_cur+win_size):
110 if not self.density_dic[chr].has_key(k_cur):
111 self.density_dic[chr][k_cur] = 0
112 self.density_dic[chr][k_cur] += score
116 k_prev = tmp_k_cur - step
117 if not chr_beg in range(k_prev, k_prev+win_size):
119 if not self.density_dic[chr].has_key(k_prev):
120 self.density_dic[chr][k_prev] = 0
121 self.density_dic[chr][k_prev] += score
126 k_next = tmp_k_cur + step
127 if not chr_beg in range(k_next, k_next+win_size):
130 if not self.density_dic[chr].has_key(k_next):
131 self.density_dic[chr][k_next] = 0
132 self.density_dic[chr][k_next] += score
136 ###########################################
138 ###########################################
139 def put_record(self, r_num):
140 rec = self.rec_dic[r_num]
142 self.out_stream.write("%s: %s\n" % (k, rec[k]))
144 ###########################################
146 ###########################################
147 def put_new_record(self):
148 chrs = self.density_dic.keys()
151 max_k = max(self.density_dic[chr].keys())
152 for k in range(0, max_k+self.step, step):
153 chr_beg, chr_end = k, k+self.win_size
155 count = self.density_dic[chr][k]
158 self.out_stream.write("WIN_CHR: %s\n" % (chr))
159 self.out_stream.write("WIN_CHR_BEG: %d\n" % (chr_beg))
160 self.out_stream.write("WIN_CHR_END: %d\n" % (chr_end))
161 self.out_stream.write("WIN_DENSITY: %d\n" % (count))
162 self.out_stream.write("---\n")
163 ###########################################
165 ###########################################
166 def print_usage(self, opt):
168 bp_dir = os.environ['BP_DIR']
169 usage_path = bp_dir + os.path.sep + "bp_usage" + os.path.sep + "calc_density.wiki"
170 os.system("print_wiki -i %s %s" % (usage_path, opt))
171 ###########################################
179 opts, args = getopt.getopt(sys.argv[1:], "I:O:?vxw:s:", ["stream_in=", "stream_out=", "help", "verbose", "no_stream", "window_size=", "step_size="])
180 except getopt.GetoptError, err:
181 # print help information and exit:
182 print str(err) # will print something like "option -a not recognized"
187 if os.isatty( sys.stdin.fileno() ):
199 if o in ("-I", "--stream_in"):
201 elif o in ("-O", "--stream_out"):
204 den.print_usage("-?")
207 den.print_usage("-?")
209 elif o in ("-v", "--verbose"):
211 elif o in ("-w", "--window_size"):
212 win_size = string.atoi(a)
214 elif o in ("-s", "--step_size"):
215 step = string.atoi(a)
219 if not win_size_flag:
222 step = win_size # default increment step is win_size
224 if not den.set_values(win_size, step):
225 sys.stderr.write("ERROR!\n")
226 sys.stderr.write("increment step must not be larger than window size\n")
230 den.open_streams(stream_in, stream_out)
232 sys.stderr.write("%s\n" % ("IOError"))
236 if not den.in_stream:
242 rec = den.get_record()
245 rec_num = den.process_record(rec)
246 den.put_record(rec_num)
247 den.put_new_record() # print out newly generated records. i.e., WIN_CHR, WIN_CHR_BEG, WIN_CHR_END and WIN_DENSITY for this script.