Mercurial > repos > xuebing > sharplabtool

comparison tools/samtools/bam_to_sam.py @ 0:9071e359b9a3

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