Mercurial > repos > iuc > genetrack

comparison genetrack_util.py @ 1:df7ac50ade5d draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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')