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