Mercurial > repos > mvdbeek > dapars

comparison dapars.py @ 22:4309798bfb61 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/mvdbeek/dapars commit b1b007c561ea6c9db145c88b6b128d66ecd05e24-dirty
author mvdbeek
date Fri, 30 Oct 2015 12:30:22 -0400
parents eace1635959a
children
comparison
equal deleted inserted replaced
21:eace1635959a 22:4309798bfb61
41 help="Start search for breakpoint n nucleotides downstream of UTR start") 41 help="Start search for breakpoint n nucleotides downstream of UTR start")
42 parser.add_argument("-ct", "--coverage_threshold", required=False, type=float, default=20, 42 parser.add_argument("-ct", "--coverage_threshold", required=False, type=float, default=20,
43 help="minimum coverage in each aligment to be considered for determining breakpoints") 43 help="minimum coverage in each aligment to be considered for determining breakpoints")
44 parser.add_argument("-b", "--breakpoint_bed", required=False, type=argparse.FileType('w'), 44 parser.add_argument("-b", "--breakpoint_bed", required=False, type=argparse.FileType('w'),
45 help="Write bedfile with coordinates of breakpoint positions to supplied path.") 45 help="Write bedfile with coordinates of breakpoint positions to supplied path.")
46 parser.add_argument("-v", "--version", action='version', version='%(prog)s 0.2.3') 46 parser.add_argument("-v", "--version", action='version', version='%(prog)s 0.2.4')
47 parser.add_argument("-p", "--plot_path", default=None, required=False, type=directory_path, 47 parser.add_argument("-p", "--plot_path", default=None, required=False, type=directory_path,
48 help="If plot_path is specified will write a coverage plot for every UTR in that directory.") 48 help="If plot_path is specified will write a coverage plot for every UTR in that directory.")
49 parser.add_argument("-html", "--html_file", default=None, required=False, type=argparse.FileType('w'), 49 parser.add_argument("-html", "--html_file", default=None, required=False, type=argparse.FileType('w'),
50 help="Write an html file to the specified location. Only to be used within a galaxy wrapper") 50 help="Write an html file to the specified location. Only to be used within a galaxy wrapper")
51 return parser.parse_args() 51 return parser.parse_args()
193 "result_d":result_d} 193 "result_d":result_d}
194 pool = Pool(self.n_cpus) 194 pool = Pool(self.n_cpus)
195 tasks = [ (self.utr_coverages[utr], self.plot_path, utr, utr_d, self.coverage_weights, len(self.control_alignments), 195 tasks = [ (self.utr_coverages[utr], self.plot_path, utr, utr_d, self.coverage_weights, len(self.control_alignments),
196 len(self.treatment_alignments), self.search_start, self.local_minimum, self.coverage_threshold) \ 196 len(self.treatment_alignments), self.search_start, self.local_minimum, self.coverage_threshold) \
197 for utr, utr_d in self.utr_dict.iteritems() ] 197 for utr, utr_d in self.utr_dict.iteritems() ]
198 #processed_tasks = [ pool.apply_async(calculate_all_utr, t) for t in tasks] 198 processed_tasks = [ pool.apply_async(calculate_all_utr, t) for t in tasks]
199 #result_list = [res.get() for res in processed_tasks] 199 result_list = [res.get() for res in processed_tasks]
200 result_list = [calculate_all_utr(*t) for t in tasks] # uncomment for easier debugging 200 #result_list = [calculate_all_utr(*t) for t in tasks] # uncomment for easier debugging
201 for res_control, res_treatment in result_list: 201 for res_control, res_treatment in result_list:
202 if not res_control: 202 if not res_control:
203 continue 203 continue
204 for i, result in enumerate(res_control): 204 for i, result in enumerate(res_control):
205 if isinstance(result, dict): 205 if isinstance(result, dict):
352 """ 352 """
353 mse_control = np.array([mse[0] for mse in mse_list]) 353 mse_control = np.array([mse[0] for mse in mse_list])
354 mse_treatment = np.array([mse[1] for mse in mse_list]) 354 mse_treatment = np.array([mse[1] for mse in mse_list])
355 control_breakpoints = list(get_minima(mse_control, local_minimum)) 355 control_breakpoints = list(get_minima(mse_control, local_minimum))
356 treatment_breakpoints = list(get_minima(mse_treatment, local_minimum)) 356 treatment_breakpoints = list(get_minima(mse_treatment, local_minimum))
357 print control_breakpoints
358 print treatment_breakpoints
359 control_abundances = [estimate_abundance(normalized_utr_coverage, bp, num_samples) for bp in control_breakpoints] 357 control_abundances = [estimate_abundance(normalized_utr_coverage, bp, num_samples) for bp in control_breakpoints]
360 treatment_abundances = [estimate_abundance(normalized_utr_coverage, bp, num_samples) for bp in treatment_breakpoints] 358 treatment_abundances = [estimate_abundance(normalized_utr_coverage, bp, num_samples) for bp in treatment_breakpoints]
361 return control_breakpoints, control_abundances, treatment_breakpoints, treatment_abundances 359 return control_breakpoints, control_abundances, treatment_breakpoints, treatment_abundances
362 360
363 def get_minima(a, local_minimum=False): 361 def get_minima(a, local_minimum=False):
364 """ 362 """
365 get minima for numpy array a. If local is false, only return absolute minimum 363 get minima for numpy array a. If local is false, only return absolute minimum
366 """ 364 """
367 if not local_minimum: 365 if not local_minimum:
368 print np.where(a == a.min())[0]
369 return np.where(a == a.min())[0] 366 return np.where(a == a.min())[0]
370 else: 367 else:
371 return np.where(np.r_[True, a[1:] < a[:-1]] & np.r_[a[:-1] < a[1:], True])[0]+1 368 return np.where(np.r_[True, a[1:] < a[:-1]] & np.r_[a[:-1] < a[1:], True])[0]+1
372 369
373 def estimate_mse(cov, bp, num_samples, num_control): 370 def estimate_mse(cov, bp, num_samples, num_control):