/tools/samtools/sam_to_bam.py

https://bitbucket.org/cistrome/cistrome-harvard/ · Python · 197 lines · 150 code · 7 blank · 40 comment · 49 complexity · 6b8395f02f8d20ed800e18b6c452dbb2 MD5 · raw file

  1. #!/usr/bin/env python
  2. """
  3. Converts SAM data to sorted BAM data.
  4. usage: sam_to_bam.py [options]
  5. --input1: SAM file to be converted
  6. --dbkey: dbkey value
  7. --ref_file: Reference file if choosing from history
  8. --output1: output dataset in bam format
  9. --index_dir: GALAXY_DATA_INDEX_DIR
  10. """
  11. import optparse, os, sys, subprocess, tempfile, shutil, gzip
  12. from galaxy import eggs
  13. import pkg_resources; pkg_resources.require( "bx-python" )
  14. from bx.cookbook import doc_optparse
  15. from galaxy import util
  16. def stop_err( msg ):
  17. sys.stderr.write( '%s\n' % msg )
  18. sys.exit()
  19. def check_seq_file( dbkey, cached_seqs_pointer_file ):
  20. seq_path = ''
  21. for line in open( cached_seqs_pointer_file ):
  22. line = line.rstrip( '\r\n' )
  23. if line and not line.startswith( '#' ) and line.startswith( 'index' ):
  24. fields = line.split( '\t' )
  25. if len( fields ) < 3:
  26. continue
  27. if fields[1] == dbkey:
  28. seq_path = fields[2].strip()
  29. break
  30. return seq_path
  31. def __main__():
  32. #Parse Command Line
  33. parser = optparse.OptionParser()
  34. parser.add_option( '', '--input1', dest='input1', help='The input SAM dataset' )
  35. parser.add_option( '', '--dbkey', dest='dbkey', help='The build of the reference dataset' )
  36. parser.add_option( '', '--ref_file', dest='ref_file', help='The reference dataset from the history' )
  37. parser.add_option( '', '--output1', dest='output1', help='The output BAM dataset' )
  38. parser.add_option( '', '--index_dir', dest='index_dir', help='GALAXY_DATA_INDEX_DIR' )
  39. ( options, args ) = parser.parse_args()
  40. # output version # of tool
  41. try:
  42. tmp = tempfile.NamedTemporaryFile().name
  43. tmp_stdout = open( tmp, 'wb' )
  44. proc = subprocess.Popen( args='samtools 2>&1', shell=True, stdout=tmp_stdout )
  45. tmp_stdout.close()
  46. returncode = proc.wait()
  47. stdout = None
  48. for line in open( tmp_stdout.name, 'rb' ):
  49. if line.lower().find( 'version' ) >= 0:
  50. stdout = line.strip()
  51. break
  52. if stdout:
  53. sys.stdout.write( 'Samtools %s\n' % stdout )
  54. else:
  55. raise Exception
  56. except:
  57. sys.stdout.write( 'Could not determine Samtools version\n' )
  58. cached_seqs_pointer_file = '%s/sam_fa_indices.loc' % options.index_dir
  59. if not os.path.exists( cached_seqs_pointer_file ):
  60. stop_err( 'The required file (%s) does not exist.' % cached_seqs_pointer_file )
  61. # If found for the dbkey, seq_path will look something like /galaxy/data/equCab2/sam_index/equCab2.fa,
  62. # and the equCab2.fa file will contain fasta sequences.
  63. seq_path = check_seq_file( options.dbkey, cached_seqs_pointer_file )
  64. tmp_dir = tempfile.mkdtemp()
  65. if not options.ref_file or options.ref_file == 'None':
  66. # We're using locally cached reference sequences( e.g., /galaxy/data/equCab2/sam_index/equCab2.fa ).
  67. # The indexes for /galaxy/data/equCab2/sam_index/equCab2.fa will be contained in
  68. # a file named /galaxy/data/equCab2/sam_index/equCab2.fa.fai
  69. fai_index_file_base = seq_path
  70. fai_index_file_path = '%s.fai' % seq_path
  71. if not os.path.exists( fai_index_file_path ):
  72. #clean up temp files
  73. if os.path.exists( tmp_dir ):
  74. shutil.rmtree( tmp_dir )
  75. stop_err( 'No sequences are available for build (%s), request them by reporting this error.' % options.dbkey )
  76. else:
  77. try:
  78. # Create indexes for history reference ( e.g., ~/database/files/000/dataset_1.dat ) using samtools faidx, which will:
  79. # - index reference sequence in the FASTA format or extract subsequence from indexed reference sequence
  80. # - if no region is specified, faidx will index the file and create <ref.fasta>.fai on the disk
  81. # - if regions are specified, the subsequences will be retrieved and printed to stdout in the FASTA format
  82. # - the input file can be compressed in the RAZF format.
  83. # IMPORTANT NOTE: a real weakness here is that we are creating indexes for the history dataset
  84. # every time we run this tool. It would be nice if we could somehow keep track of user's specific
  85. # index files so they could be re-used.
  86. fai_index_file_base = tempfile.NamedTemporaryFile( dir=tmp_dir ).name
  87. # At this point, fai_index_file_path will look something like /tmp/dataset_13.dat
  88. os.symlink( options.ref_file, fai_index_file_base )
  89. fai_index_file_path = '%s.fai' % fai_index_file_base
  90. command = 'samtools faidx %s' % fai_index_file_base
  91. tmp = tempfile.NamedTemporaryFile( dir=tmp_dir ).name
  92. tmp_stderr = open( tmp, 'wb' )
  93. proc = subprocess.Popen( args=command, shell=True, cwd=tmp_dir, stderr=tmp_stderr.fileno() )
  94. returncode = proc.wait()
  95. tmp_stderr.close()
  96. # get stderr, allowing for case where it's very large
  97. tmp_stderr = open( tmp, 'rb' )
  98. stderr = ''
  99. buffsize = 1048576
  100. try:
  101. while True:
  102. stderr += tmp_stderr.read( buffsize )
  103. if not stderr or len( stderr ) % buffsize != 0:
  104. break
  105. except OverflowError:
  106. pass
  107. tmp_stderr.close()
  108. if returncode != 0:
  109. raise Exception, stderr
  110. if os.path.getsize( fai_index_file_path ) == 0:
  111. raise Exception, 'Index file empty, there may be an error with your reference file or settings.'
  112. except Exception, e:
  113. #clean up temp files
  114. if os.path.exists( tmp_dir ):
  115. shutil.rmtree( tmp_dir )
  116. stop_err( 'Error creating indexes from reference (%s), %s' % ( options.ref_file, str( e ) ) )
  117. try:
  118. # Extract all alignments from the input SAM file to BAM format ( since no region is specified, all the alignments will be extracted ).
  119. tmp_aligns_file = tempfile.NamedTemporaryFile( dir=tmp_dir )
  120. tmp_aligns_file_name = tmp_aligns_file.name
  121. tmp_aligns_file.close()
  122. command = 'samtools view -bt %s -o %s %s' % ( fai_index_file_path, tmp_aligns_file_name, options.input1 )
  123. tmp = tempfile.NamedTemporaryFile( dir=tmp_dir ).name
  124. tmp_stderr = open( tmp, 'wb' )
  125. proc = subprocess.Popen( args=command, shell=True, cwd=tmp_dir, stderr=tmp_stderr.fileno() )
  126. returncode = proc.wait()
  127. tmp_stderr.close()
  128. # get stderr, allowing for case where it's very large
  129. tmp_stderr = open( tmp, 'rb' )
  130. stderr = ''
  131. buffsize = 1048576
  132. try:
  133. while True:
  134. stderr += tmp_stderr.read( buffsize )
  135. if not stderr or len( stderr ) % buffsize != 0:
  136. break
  137. except OverflowError:
  138. pass
  139. tmp_stderr.close()
  140. if returncode != 0:
  141. raise Exception, stderr
  142. except Exception, e:
  143. #clean up temp files
  144. if os.path.exists( tmp_dir ):
  145. shutil.rmtree( tmp_dir )
  146. stop_err( 'Error extracting alignments from (%s), %s' % ( options.input1, str( e ) ) )
  147. try:
  148. # Sort alignments by leftmost coordinates. File <out.prefix>.bam will be created. This command
  149. # may also create temporary files <out.prefix>.%d.bam when the whole alignment cannot be fitted
  150. # into memory ( controlled by option -m ).
  151. tmp_sorted_aligns_file = tempfile.NamedTemporaryFile( dir=tmp_dir )
  152. tmp_sorted_aligns_file_name = tmp_sorted_aligns_file.name
  153. tmp_sorted_aligns_file.close()
  154. command = 'samtools sort %s %s' % ( tmp_aligns_file_name, tmp_sorted_aligns_file_name )
  155. tmp = tempfile.NamedTemporaryFile( dir=tmp_dir ).name
  156. tmp_stderr = open( tmp, 'wb' )
  157. proc = subprocess.Popen( args=command, shell=True, cwd=tmp_dir, stderr=tmp_stderr.fileno() )
  158. returncode = proc.wait()
  159. tmp_stderr.close()
  160. # get stderr, allowing for case where it's very large
  161. tmp_stderr = open( tmp, 'rb' )
  162. stderr = ''
  163. buffsize = 1048576
  164. try:
  165. while True:
  166. stderr += tmp_stderr.read( buffsize )
  167. if not stderr or len( stderr ) % buffsize != 0:
  168. break
  169. except OverflowError:
  170. pass
  171. tmp_stderr.close()
  172. if returncode != 0:
  173. raise Exception, stderr
  174. except Exception, e:
  175. #clean up temp files
  176. if os.path.exists( tmp_dir ):
  177. shutil.rmtree( tmp_dir )
  178. stop_err( 'Error sorting alignments from (%s), %s' % ( tmp_aligns_file_name, str( e ) ) )
  179. # Move tmp_aligns_file_name to our output dataset location
  180. sorted_bam_file = '%s.bam' % tmp_sorted_aligns_file_name
  181. shutil.move( sorted_bam_file, options.output1 )
  182. #clean up temp files
  183. if os.path.exists( tmp_dir ):
  184. shutil.rmtree( tmp_dir )
  185. # check that there are results in the output file
  186. if os.path.getsize( options.output1 ) > 0:
  187. sys.stdout.write( 'SAM file converted to BAM' )
  188. else:
  189. stop_err( 'Error creating sorted version of BAM file.' )
  190. if __name__=="__main__": __main__()