Mercurial > repos > iuc > genetrack
comparison genetrack_util.py @ 1:df7ac50ade5d draft
planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/genetrack commit 6cd0f93f4eb8649802683ae4c189c8ad48827a49
author | iuc |
---|---|
date | Fri, 13 Jan 2017 10:20:48 -0500 |
parents | 25cd59a002d9 |
children | 41887967ef14 |
comparison
equal
deleted
inserted
replaced
0:25cd59a002d9 | 1:df7ac50ade5d |
---|---|
255 should be at). | 255 should be at). |
256 """ | 256 """ |
257 lo, hi = get_range(data) | 257 lo, hi = get_range(data) |
258 rng = hi - lo | 258 rng = hi - lo |
259 shift = width - lo | 259 shift = width - lo |
260 return numpy.zeros(rng+width*2, numpy.float), shift | 260 return numpy.zeros(rng + width * 2, numpy.float), shift |
261 | 261 |
262 | 262 |
263 def normal_array(width, sigma, normalize=True): | 263 def normal_array(width, sigma, normalize=True): |
264 """ | 264 """ |
265 Returns an array of the normal distribution of the specified width | 265 Returns an array of the normal distribution of the specified width |
272 # width is the half of the distribution | 272 # width is the half of the distribution |
273 values = map(normal_func, range(-width, width)) | 273 values = map(normal_func, range(-width, width)) |
274 values = numpy.array(values, numpy.float) | 274 values = numpy.array(values, numpy.float) |
275 # normalization | 275 # normalization |
276 if normalize: | 276 if normalize: |
277 values = 1.0/math.sqrt(2 * numpy.pi * sigma2) * values | 277 values = 1.0 / math.sqrt(2 * numpy.pi * sigma2) * values |
278 return values | 278 return values |
279 | 279 |
280 | 280 |
281 def call_peaks(array, shift, data, keys, direction, down_width, up_width, exclusion): | 281 def call_peaks(array, shift, data, keys, direction, down_width, up_width, exclusion): |
282 peaks = [] | 282 peaks = [] |
290 neg = up_width or exclusion // 2 | 290 neg = up_width or exclusion // 2 |
291 # Reverse strand | 291 # Reverse strand |
292 if direction == 2: | 292 if direction == 2: |
293 # Swap positive and negative widths | 293 # Swap positive and negative widths |
294 pos, neg = neg, pos | 294 pos, neg = neg, pos |
295 peaks.append(Peak(int(index)-shift, pos, neg)) | 295 peaks.append(Peak(int(index) - shift, pos, neg)) |
296 find_peaks() | 296 find_peaks() |
297 | 297 |
298 def calculate_reads(): | 298 def calculate_reads(): |
299 # Calculate the number of reads in each peak | 299 # Calculate the number of reads in each peak |
300 for peak in peaks: | 300 for peak in peaks: |
311 peaks_by_value = peaks[:] | 311 peaks_by_value = peaks[:] |
312 peaks_by_value.sort(key=lambda peak: -peak.value) | 312 peaks_by_value.sort(key=lambda peak: -peak.value) |
313 for peak in peaks_by_value: | 313 for peak in peaks_by_value: |
314 peak.safe = True | 314 peak.safe = True |
315 window = get_window(peaks, | 315 window = get_window(peaks, |
316 peak.index-exclusion//2, | 316 peak.index - exclusion // 2, |
317 peak.index+exclusion//2, | 317 peak.index + exclusion // 2, |
318 peak_keys) | 318 peak_keys) |
319 for excluded in window: | 319 for excluded in window: |
320 if excluded.safe: | 320 if excluded.safe: |
321 continue | 321 continue |
322 i = get_index(excluded.index, peak_keys) | 322 i = get_index(excluded.index, peak_keys) |
344 # Add each read's normal to the array | 344 # Add each read's normal to the array |
345 for read in data: | 345 for read in data: |
346 index, forward, reverse = read | 346 index, forward, reverse = read |
347 # Add the normals to the appropriate regions | 347 # Add the normals to the appropriate regions |
348 if forward: | 348 if forward: |
349 forward_array[index+forward_shift-width:index+forward_shift+width] += normal * forward | 349 forward_array[index + forward_shift - width:index + forward_shift + width] += normal * forward |
350 if reverse: | 350 if reverse: |
351 reverse_array[index+reverse_shift-width:index+reverse_shift+width] += normal * reverse | 351 reverse_array[index + reverse_shift - width:index + reverse_shift + width] += normal * reverse |
352 populate_array() | 352 populate_array() |
353 forward_peaks = call_peaks(forward_array, forward_shift, data, keys, 1, down_width, up_width, exclusion) | 353 forward_peaks = call_peaks(forward_array, forward_shift, data, keys, 1, down_width, up_width, exclusion) |
354 reverse_peaks = call_peaks(reverse_array, reverse_shift, data, keys, 2, down_width, up_width, exclusion) | 354 reverse_peaks = call_peaks(reverse_array, reverse_shift, data, keys, 2, down_width, up_width, exclusion) |
355 # Convert chromosome name in preparation for writing output | 355 # Convert chromosome name in preparation for writing output |
356 cname = convert_data(cname, 'zeropad', 'numeric') | 356 cname = convert_data(cname, 'zeropad', 'numeric') |