Mercurial > repos > bgruening > bismark

comparison bismark_wrapper.py @ 4:243e8f9fb75b draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author bgruening
date Mon, 09 Feb 2015 18:24:41 -0500
parents 91f07ff056ca
children
comparison
equal deleted inserted replaced
3:91f07ff056ca 4:243e8f9fb75b
15 sys.stderr.write( "%s\n" % msg ) 15 sys.stderr.write( "%s\n" % msg )
16 sys.exit() 16 sys.exit()
17 17
18 def __main__(): 18 def __main__():
19 19
20 print 'tempfile_location',tempfile.gettempdir()
21 #Parse Command Line 20 #Parse Command Line
22 parser = argparse.ArgumentParser(description='Wrapper for the bismark bisulfite mapper.') 21 parser = argparse.ArgumentParser(description='Wrapper for the bismark bisulfite mapper.')
23 parser.add_argument( '-p', '--num-threads', dest='num_threads', 22 parser.add_argument( '-p', '--num-threads', dest='num_threads',
24 type=int, default=4, help='Use this many threads to align reads. The default is 4.' ) 23 type=int, default=4, help='Use this many threads to align reads. The default is 4.' )
25 24
256 255
257 arguments.update( {'args': additional_opts, 'reads': reads} ) 256 arguments.update( {'args': additional_opts, 'reads': reads} )
258 257
259 # Final bismark command: 258 # Final bismark command:
260 cmd = cmd % arguments 259 cmd = cmd % arguments
261 print 'bismark_cmd:', cmd
262 #sys.stderr.write( cmd )
263 #sys.exit(1)
264 # Run 260 # Run
265 try: 261 try:
266 tmp_out = tempfile.NamedTemporaryFile().name 262 tmp_out = tempfile.NamedTemporaryFile().name
267 tmp_stdout = open( tmp_out, 'wb' ) 263 tmp_stdout = open( tmp_out, 'wb' )
268 tmp_err = tempfile.NamedTemporaryFile().name 264 tmp_err = tempfile.NamedTemporaryFile().name
319 """ 315 """
320 merge all bam files 316 merge all bam files
321 """ 317 """
322 #tmp_out = tempfile.NamedTemporaryFile( dir=output_dir ).name 318 #tmp_out = tempfile.NamedTemporaryFile( dir=output_dir ).name
323 tmp_stdout = open( tmp_out, 'wab' ) 319 tmp_stdout = open( tmp_out, 'wab' )
324 tmp_err = tempfile.NamedTemporaryFile( dir=output_dir ).name 320 #tmp_err = tempfile.NamedTemporaryFile( dir=output_dir ).name
325 tmp_stderr = open( tmp_err, 'wb' ) 321 tmp_stderr = open( tmp_err, 'wab' )
326 322
327 tmp_res = tempfile.NamedTemporaryFile( dir= output_dir).name 323 tmp_res = tempfile.NamedTemporaryFile( dir= output_dir).name
328 324
329 bam_files = glob( os.path.join( output_dir, '*.bam') ) 325 bam_files = glob( os.path.join( output_dir, '*.bam') )
330 if len( bam_files ) > 1: 326 if len( bam_files ) > 1:
336 tmp_stdout.close() 332 tmp_stdout.close()
337 tmp_stderr.close() 333 tmp_stderr.close()
338 if returncode != 0: 334 if returncode != 0:
339 raise Exception, open( tmp_stderr.name ).read() 335 raise Exception, open( tmp_stderr.name ).read()
340 else: 336 else:
341 tmp_res = bam_files[0] 337 tmp_res = bam_files[0]
342 338
343 bam_path = "%s" % tmp_res 339 bam_path = "%s" % tmp_res
344 340
345 if os.path.exists( bam_path ): 341 if os.path.exists( bam_path ):
346 if args.sort_bam: 342 if args.sort_bam:
347 cmd = 'samtools sort -@ %s %s %s' % (args.num_threads, bam_path, args.output) 343 cmd = 'samtools sort -@ %s %s sorted_bam' % (args.num_threads, bam_path)
348 else: 344 proc = subprocess.Popen( args=shlex.split( cmd ) )
349 shutil.copy( bam_path, args.output ) 345 returncode = proc.wait()
346 if returncode != 0:
347 raise Exception("Error during '%s'" % cmd)
348 shutil.move( 'sorted_bam.bam', args.output )
349 else:
350 shutil.move( bam_path, args.output )
350 else: 351 else:
351 stop_err( 'BAM file no found:\n' + str( bam_path ) ) 352 stop_err( 'BAM file no found:\n' + str( bam_path ) )
352
353 353
354 354
355 # TODO: look for errors in program output. 355 # TODO: look for errors in program output.
356 except Exception, e: 356 except Exception, e:
357 stop_err( 'Error in merging bam files:\n' + str( e ) ) 357 stop_err( 'Error in merging bam files:\n' + str( e ) )