Mercurial > repos > devteam > variant_filtration

comparison gatk_wrapper.py @ 0:da6e2503c62d draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Imported from capsule None
author devteam
date Tue, 01 Apr 2014 10:50:10 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:da6e2503c62d
1 #!/usr/bin/env python
2 #Dan Blankenberg
3
4 """
5 A wrapper script for running the GenomeAnalysisTK.jar commands.
6 """
7
8 import sys, optparse, os, tempfile, subprocess, shutil
9 from binascii import unhexlify
10 from string import Template
11
12 GALAXY_EXT_TO_GATK_EXT = { 'gatk_interval':'intervals', 'bam_index':'bam.bai', 'gatk_dbsnp':'dbSNP', 'picard_interval_list':'interval_list' } #items not listed here will use the galaxy extension as-is
13 GALAXY_EXT_TO_GATK_FILE_TYPE = GALAXY_EXT_TO_GATK_EXT #for now, these are the same, but could be different if needed
14 DEFAULT_GATK_PREFIX = "gatk_file"
15 CHUNK_SIZE = 2**20 #1mb
16
17
18 def cleanup_before_exit( tmp_dir ):
19 if tmp_dir and os.path.exists( tmp_dir ):
20 shutil.rmtree( tmp_dir )
21
22 def gatk_filename_from_galaxy( galaxy_filename, galaxy_ext, target_dir = None, prefix = None ):
23 suffix = GALAXY_EXT_TO_GATK_EXT.get( galaxy_ext, galaxy_ext )
24 if prefix is None:
25 prefix = DEFAULT_GATK_PREFIX
26 if target_dir is None:
27 target_dir = os.getcwd()
28 gatk_filename = os.path.join( target_dir, "%s.%s" % ( prefix, suffix ) )
29 os.symlink( galaxy_filename, gatk_filename )
30 return gatk_filename
31
32 def gatk_filetype_argument_substitution( argument, galaxy_ext ):
33 return argument % dict( file_type = GALAXY_EXT_TO_GATK_FILE_TYPE.get( galaxy_ext, galaxy_ext ) )
34
35 def open_file_from_option( filename, mode = 'rb' ):
36 if filename:
37 return open( filename, mode = mode )
38 return None
39
40 def html_report_from_directory( html_out, dir ):
41 html_out.write( '<html>\n<head>\n<title>Galaxy - GATK Output</title>\n</head>\n<body>\n<p/>\n<ul>\n' )
42 for fname in sorted( os.listdir( dir ) ):
43 html_out.write( '<li><a href="%s">%s</a></li>\n' % ( fname, fname ) )
44 html_out.write( '</ul>\n</body>\n</html>\n' )
45
46 def index_bam_files( bam_filenames, tmp_dir ):
47 for bam_filename in bam_filenames:
48 bam_index_filename = "%s.bai" % bam_filename
49 if not os.path.exists( bam_index_filename ):
50 #need to index this bam file
51 stderr_name = tempfile.NamedTemporaryFile( prefix = "bam_index_stderr" ).name
52 command = 'samtools index %s %s' % ( bam_filename, bam_index_filename )
53 proc = subprocess.Popen( args=command, shell=True, stderr=open( stderr_name, 'wb' ) )
54 return_code = proc.wait()
55 if return_code:
56 for line in open( stderr_name ):
57 print >> sys.stderr, line
58 os.unlink( stderr_name ) #clean up
59 cleanup_before_exit( tmp_dir )
60 raise Exception( "Error indexing BAM file" )
61 os.unlink( stderr_name ) #clean up
62
63 def __main__():
64 #Parse Command Line
65 parser = optparse.OptionParser()
66 parser.add_option( '-p', '--pass_through', dest='pass_through_options', action='append', type="string", help='These options are passed through directly to GATK, without any modification.' )
67 parser.add_option( '-o', '--pass_through_options', dest='pass_through_options_encoded', action='append', type="string", help='These options are passed through directly to GATK, with decoding from binascii.unhexlify.' )
68 parser.add_option( '-d', '--dataset', dest='datasets', action='append', type="string", nargs=4, help='"-argument" "original_filename" "galaxy_filetype" "name_prefix"' )
69 parser.add_option( '', '--max_jvm_heap', dest='max_jvm_heap', action='store', type="string", default=None, help='If specified, the maximum java virtual machine heap size will be set to the provide value.' )
70 parser.add_option( '', '--max_jvm_heap_fraction', dest='max_jvm_heap_fraction', action='store', type="int", default=None, help='If specified, the maximum java virtual machine heap size will be set to the provide value as a fraction of total physical memory.' )
71 parser.add_option( '', '--stdout', dest='stdout', action='store', type="string", default=None, help='If specified, the output of stdout will be written to this file.' )
72 parser.add_option( '', '--stderr', dest='stderr', action='store', type="string", default=None, help='If specified, the output of stderr will be written to this file.' )
73 parser.add_option( '', '--html_report_from_directory', dest='html_report_from_directory', action='append', type="string", nargs=2, help='"Target HTML File" "Directory"')
74 (options, args) = parser.parse_args()
75
76 tmp_dir = tempfile.mkdtemp( prefix='tmp-gatk-' )
77 if options.pass_through_options:
78 cmd = ' '.join( options.pass_through_options )
79 else:
80 cmd = ''
81 if options.pass_through_options_encoded:
82 cmd = '%s %s' % ( cmd, ' '.join( map( unhexlify, options.pass_through_options_encoded ) ) )
83 if options.max_jvm_heap is not None:
84 cmd = cmd.replace( 'java ', 'java -Xmx%s ' % ( options.max_jvm_heap ), 1 )
85 elif options.max_jvm_heap_fraction is not None:
86 cmd = cmd.replace( 'java ', 'java -XX:DefaultMaxRAMFraction=%s -XX:+UseParallelGC ' % ( options.max_jvm_heap_fraction ), 1 )
87 bam_filenames = []
88 if options.datasets:
89 for ( dataset_arg, filename, galaxy_ext, prefix ) in options.datasets:
90 gatk_filename = gatk_filename_from_galaxy( filename, galaxy_ext, target_dir = tmp_dir, prefix = prefix )
91 if dataset_arg:
92 cmd = '%s %s "%s"' % ( cmd, gatk_filetype_argument_substitution( dataset_arg, galaxy_ext ), gatk_filename )
93 if galaxy_ext == "bam":
94 bam_filenames.append( gatk_filename )
95 index_bam_files( bam_filenames, tmp_dir )
96 #set up stdout and stderr output options
97 stdout = open_file_from_option( options.stdout, mode = 'wb' )
98 stderr = open_file_from_option( options.stderr, mode = 'wb' )
99 #if no stderr file is specified, we'll use our own
100 if stderr is None:
101 stderr = tempfile.NamedTemporaryFile( prefix="gatk-stderr-", dir=tmp_dir )
102
103 proc = subprocess.Popen( args=cmd, stdout=stdout, stderr=stderr, shell=True, cwd=tmp_dir )
104 return_code = proc.wait()
105
106 if return_code:
107 stderr_target = sys.stderr
108 else:
109 stderr_target = sys.stdout
110 stderr.flush()
111 stderr.seek(0)
112 while True:
113 chunk = stderr.read( CHUNK_SIZE )
114 if chunk:
115 stderr_target.write( chunk )
116 else:
117 break
118 stderr.close()
119 #generate html reports
120 if options.html_report_from_directory:
121 for ( html_filename, html_dir ) in options.html_report_from_directory:
122 html_report_from_directory( open( html_filename, 'wb' ), html_dir )
123
124 cleanup_before_exit( tmp_dir )
125
126 if __name__=="__main__": __main__()