/tools/gatk/gatk_wrapper.py

https://bitbucket.org/cistrome/cistrome-harvard/ · Python · 126 lines · 100 code · 16 blank · 10 comment · 28 complexity · 90f5a44f1e756b8e3020fc40fb22ed76 MD5 · raw file 

    

  1. #!/usr/bin/env python
    2. 

  2. #Dan Blankenberg
    3. 

  3. 

  4. """
    5. 

  5. A wrapper script for running the GenomeAnalysisTK.jar commands.
    6. 

  6. """
    7. 

    8. 

  8. import sys, optparse, os, tempfile, subprocess, shutil
    9. 

  9. from binascii import unhexlify
    10. 

  10. from string import Template
    11. 

  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. 

  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. 

  14. DEFAULT_GATK_PREFIX = "gatk_file"
    15. 

  15. CHUNK_SIZE = 2**20 #1mb
    16. 

  16. 

  17. 

  18. def cleanup_before_exit( tmp_dir ):
    19. 

  19.     if tmp_dir and os.path.exists( tmp_dir ):
    20. 

  20.         shutil.rmtree( tmp_dir )
    21. 

  21. 

  22. def gatk_filename_from_galaxy( galaxy_filename, galaxy_ext, target_dir = None, prefix = None ):
    23. 

  23.     suffix = GALAXY_EXT_TO_GATK_EXT.get( galaxy_ext, galaxy_ext )
    24. 

  24.     if prefix is None:
    25. 

  25.         prefix = DEFAULT_GATK_PREFIX
    26. 

  26.     if target_dir is None:
    27. 

  27.         target_dir = os.getcwd()
    28. 

  28.     gatk_filename = os.path.join( target_dir, "%s.%s" % ( prefix, suffix ) )
    29. 

  29.     os.symlink( galaxy_filename, gatk_filename )
    30. 

  30.     return gatk_filename
    31. 

  31. 

  32. def gatk_filetype_argument_substitution( argument, galaxy_ext ):
    33. 

  33.     return argument % dict( file_type = GALAXY_EXT_TO_GATK_FILE_TYPE.get( galaxy_ext, galaxy_ext ) )
    34. 

  34. 

  35. def open_file_from_option( filename, mode = 'rb' ):
    36. 

  36.     if filename:
    37. 

  37.         return open( filename, mode = mode )
    38. 

  38.     return None
    39. 

  39. 

  40. def html_report_from_directory( html_out, dir ):
    41. 

  41.     html_out.write( '<html>\n<head>\n<title>Galaxy - GATK Output</title>\n</head>\n<body>\n<p/>\n<ul>\n' )
    42. 

  42.     for fname in sorted( os.listdir( dir ) ):
    43. 

  43.         html_out.write(  '<li><a href="%s">%s</a></li>\n' % ( fname, fname ) )
    44. 

  44.     html_out.write( '</ul>\n</body>\n</html>\n' )
    45. 

  45. 

  46. def index_bam_files( bam_filenames, tmp_dir ):
    47. 

  47.     for bam_filename in bam_filenames:
    48. 

  48.         bam_index_filename = "%s.bai" % bam_filename
    49. 

  49.         if not os.path.exists( bam_index_filename ):
    50. 

  50.             #need to index this bam file
    51. 

  51.             stderr_name = tempfile.NamedTemporaryFile( prefix = "bam_index_stderr" ).name
    52. 

  52.             command = 'samtools index %s %s' % ( bam_filename, bam_index_filename )
    53. 

  53.             proc = subprocess.Popen( args=command, shell=True, stderr=open( stderr_name, 'wb' ) )
    54. 

  54.             return_code = proc.wait()
    55. 

  55.             if return_code:
    56. 

  56.                 for line in open( stderr_name ):
    57. 

  57.                     print >> sys.stderr, line
    58. 

  58.                 os.unlink( stderr_name ) #clean up
    59. 

  59.                 cleanup_before_exit( tmp_dir )
    60. 

  60.                 raise Exception( "Error indexing BAM file" )
    61. 

  61.             os.unlink( stderr_name ) #clean up
    62. 

  62. 

  63. def __main__():
    64. 

  64.     #Parse Command Line
    65. 

  65.     parser = optparse.OptionParser()
    66. 

  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. 

  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. 

  68.     parser.add_option( '-d', '--dataset', dest='datasets', action='append', type="string", nargs=4, help='"-argument" "original_filename" "galaxy_filetype" "name_prefix"' )
    69. 

  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. 

  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. 

  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. 

  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. 

  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. 

  74.     (options, args) = parser.parse_args()
    75. 

  75.     
    76. 

  76.     tmp_dir = tempfile.mkdtemp( prefix='tmp-gatk-' )
    77. 

  77.     if options.pass_through_options:
    78. 

  78.         cmd = ' '.join( options.pass_through_options )
    79. 

  79.     else:
    80. 

  80.         cmd = ''
    81. 

  81.     if options.pass_through_options_encoded:
    82. 

  82.         cmd = '%s %s' % ( cmd, ' '.join( map( unhexlify, options.pass_through_options_encoded ) ) )
    83. 

  83.     if options.max_jvm_heap is not None:
    84. 

  84.         cmd = cmd.replace( 'java ', 'java -Xmx%s ' % ( options.max_jvm_heap ), 1 )
    85. 

  85.     elif options.max_jvm_heap_fraction is not None:
    86. 

  86.         cmd = cmd.replace( 'java ', 'java -XX:DefaultMaxRAMFraction=%s  -XX:+UseParallelGC ' % ( options.max_jvm_heap_fraction ), 1 )
    87. 

  87.     bam_filenames = []
    88. 

  88.     if options.datasets:
    89. 

  89.         for ( dataset_arg, filename, galaxy_ext, prefix ) in options.datasets:
    90. 

  90.             gatk_filename = gatk_filename_from_galaxy( filename, galaxy_ext, target_dir = tmp_dir, prefix = prefix )
    91. 

  91.             if dataset_arg:
    92. 

  92.                 cmd = '%s %s "%s"' % ( cmd, gatk_filetype_argument_substitution( dataset_arg, galaxy_ext ), gatk_filename )
    93. 

  93.             if galaxy_ext == "bam":
    94. 

  94.                 bam_filenames.append( gatk_filename )
    95. 

  95.     index_bam_files( bam_filenames, tmp_dir )
    96. 

  96.     #set up stdout and stderr output options
    97. 

  97.     stdout = open_file_from_option( options.stdout, mode = 'wb' )
    98. 

  98.     stderr = open_file_from_option( options.stderr, mode = 'wb' )
    99. 

  99.     #if no stderr file is specified, we'll use our own
    100. 

  100.     if stderr is None:
    101. 

  101.         stderr = tempfile.NamedTemporaryFile( prefix="gatk-stderr-", dir=tmp_dir )
    102. 

  102.     
    103. 

  103.     proc = subprocess.Popen( args=cmd, stdout=stdout, stderr=stderr, shell=True, cwd=tmp_dir )
    104. 

  104.     return_code = proc.wait()
    105. 

  105.     
    106. 

  106.     if return_code:
    107. 

  107.         stderr_target = sys.stderr
    108. 

  108.     else:
    109. 

  109.         stderr_target = sys.stdout
    110. 

  110.     stderr.flush()
    111. 

  111.     stderr.seek(0)
    112. 

  112.     while True:
    113. 

  113.         chunk = stderr.read( CHUNK_SIZE )
    114. 

  114.         if chunk:
    115. 

  115.             stderr_target.write( chunk )
    116. 

  116.         else:
    117. 

  117.             break
    118. 

  118.     stderr.close()
    119. 

  119.     #generate html reports
    120. 

  120.     if options.html_report_from_directory:
    121. 

  121.         for ( html_filename, html_dir ) in options.html_report_from_directory:
    122. 

  122.             html_report_from_directory( open( html_filename, 'wb' ), html_dir )
    123. 

  123.     
    124. 

  124.     cleanup_before_exit( tmp_dir )
    125. 

  125. 

  126. if __name__=="__main__": __main__()
    127.