/tools/samtools/bam_to_sam.py

https://bitbucket.org/cistrome/cistrome-harvard/ · Python · 129 lines · 98 code · 9 blank · 22 comment · 31 complexity · 8bef67ec54e69d1a42ac689464c63890 MD5 · raw file

  1. #!/usr/bin/env python
  2. """
  3. Converts BAM data to sorted SAM data.
  4. usage: bam_to_sam.py [options]
  5. --input1: SAM file to be converted
  6. --output1: output dataset in bam format
  7. """
  8. import optparse, os, sys, subprocess, tempfile, shutil
  9. from galaxy import eggs
  10. import pkg_resources; pkg_resources.require( "bx-python" )
  11. from bx.cookbook import doc_optparse
  12. #from galaxy import util
  13. def stop_err( msg ):
  14. sys.stderr.write( '%s\n' % msg )
  15. sys.exit()
  16. def __main__():
  17. #Parse Command Line
  18. parser = optparse.OptionParser()
  19. parser.add_option( '', '--input1', dest='input1', help='The input SAM dataset' )
  20. parser.add_option( '', '--output1', dest='output1', help='The output BAM dataset' )
  21. parser.add_option( '', '--header', dest='header', action='store_true', default=False, help='Write SAM Header' )
  22. ( options, args ) = parser.parse_args()
  23. # output version # of tool
  24. try:
  25. tmp = tempfile.NamedTemporaryFile().name
  26. tmp_stdout = open( tmp, 'wb' )
  27. proc = subprocess.Popen( args='samtools 2>&1', shell=True, stdout=tmp_stdout )
  28. tmp_stdout.close()
  29. returncode = proc.wait()
  30. stdout = None
  31. for line in open( tmp_stdout.name, 'rb' ):
  32. if line.lower().find( 'version' ) >= 0:
  33. stdout = line.strip()
  34. break
  35. if stdout:
  36. sys.stdout.write( 'Samtools %s\n' % stdout )
  37. else:
  38. raise Exception
  39. except:
  40. sys.stdout.write( 'Could not determine Samtools version\n' )
  41. tmp_dir = tempfile.mkdtemp()
  42. try:
  43. # exit if input file empty
  44. if os.path.getsize( options.input1 ) == 0:
  45. raise Exception, 'Initial BAM file empty'
  46. # Sort alignments by leftmost coordinates. File <out.prefix>.bam will be created. This command
  47. # may also create temporary files <out.prefix>.%d.bam when the whole alignment cannot be fitted
  48. # into memory ( controlled by option -m ).
  49. tmp_sorted_aligns_file = tempfile.NamedTemporaryFile( dir=tmp_dir )
  50. tmp_sorted_aligns_file_base = tmp_sorted_aligns_file.name
  51. tmp_sorted_aligns_file_name = '%s.bam' % tmp_sorted_aligns_file.name
  52. tmp_sorted_aligns_file.close()
  53. command = 'samtools sort %s %s' % ( options.input1, tmp_sorted_aligns_file_base )
  54. tmp = tempfile.NamedTemporaryFile( dir=tmp_dir ).name
  55. tmp_stderr = open( tmp, 'wb' )
  56. proc = subprocess.Popen( args=command, shell=True, cwd=tmp_dir, stderr=tmp_stderr.fileno() )
  57. returncode = proc.wait()
  58. tmp_stderr.close()
  59. # get stderr, allowing for case where it's very large
  60. tmp_stderr = open( tmp, 'rb' )
  61. stderr = ''
  62. buffsize = 1048576
  63. try:
  64. while True:
  65. stderr += tmp_stderr.read( buffsize )
  66. if not stderr or len( stderr ) % buffsize != 0:
  67. break
  68. except OverflowError:
  69. pass
  70. tmp_stderr.close()
  71. if returncode != 0:
  72. raise Exception, stderr
  73. # exit if sorted BAM file empty
  74. if os.path.getsize( tmp_sorted_aligns_file_name) == 0:
  75. raise Exception, 'Intermediate sorted BAM file empty'
  76. except Exception, e:
  77. #clean up temp files
  78. if os.path.exists( tmp_dir ):
  79. shutil.rmtree( tmp_dir )
  80. stop_err( 'Error sorting alignments from (%s), %s' % ( options.input1, str( e ) ) )
  81. try:
  82. # Extract all alignments from the input BAM file to SAM format ( since no region is specified, all the alignments will be extracted ).
  83. if options.header:
  84. view_options = "-h"
  85. else:
  86. view_options = ""
  87. command = 'samtools view %s -o %s %s' % ( view_options, options.output1, tmp_sorted_aligns_file_name )
  88. tmp = tempfile.NamedTemporaryFile( dir=tmp_dir ).name
  89. tmp_stderr = open( tmp, 'wb' )
  90. proc = subprocess.Popen( args=command, shell=True, cwd=tmp_dir, stderr=tmp_stderr.fileno() )
  91. returncode = proc.wait()
  92. tmp_stderr.close()
  93. # get stderr, allowing for case where it's very large
  94. tmp_stderr = open( tmp, 'rb' )
  95. stderr = ''
  96. buffsize = 1048576
  97. try:
  98. while True:
  99. stderr += tmp_stderr.read( buffsize )
  100. if not stderr or len( stderr ) % buffsize != 0:
  101. break
  102. except OverflowError:
  103. pass
  104. tmp_stderr.close()
  105. if returncode != 0:
  106. raise Exception, stderr
  107. except Exception, e:
  108. #clean up temp files
  109. if os.path.exists( tmp_dir ):
  110. shutil.rmtree( tmp_dir )
  111. stop_err( 'Error extracting alignments from (%s), %s' % ( options.input1, str( e ) ) )
  112. #clean up temp files
  113. if os.path.exists( tmp_dir ):
  114. shutil.rmtree( tmp_dir )
  115. # check that there are results in the output file
  116. if os.path.getsize( options.output1 ) > 0:
  117. sys.stdout.write( 'BAM file converted to SAM' )
  118. else:
  119. stop_err( 'The output file is empty, there may be an error with your input file.' )
  120. if __name__=="__main__": __main__()