|
0
|
1 #!/usr/bin/env python
|
|
|
2
|
|
|
3 """
|
|
|
4 Runs BWA on single-end or paired-end data.
|
|
|
5 Produces a SAM file containing the mappings.
|
|
|
6 Works with BWA version 0.5.9.
|
|
|
7
|
|
|
8 usage: bwa_wrapper.py [options]
|
|
|
9
|
|
|
10 See below for options
|
|
|
11 """
|
|
|
12
|
|
|
13 import optparse, os, shutil, subprocess, sys, tempfile
|
|
|
14
|
|
|
15 def stop_err( msg ):
|
|
|
16 sys.stderr.write( '%s\n' % msg )
|
|
|
17 sys.exit()
|
|
|
18
|
|
|
19 def check_is_double_encoded( fastq ):
|
|
|
20 # check that first read is bases, not one base followed by numbers
|
|
|
21 bases = [ 'A', 'C', 'G', 'T', 'a', 'c', 'g', 't', 'N' ]
|
|
|
22 nums = [ '0', '1', '2', '3' ]
|
|
|
23 for line in file( fastq, 'rb'):
|
|
|
24 if not line.strip() or line.startswith( '@' ):
|
|
|
25 continue
|
|
|
26 if len( [ b for b in line.strip() if b in nums ] ) > 0:
|
|
|
27 return False
|
|
|
28 elif line.strip()[0] in bases and len( [ b for b in line.strip() if b in bases ] ) == len( line.strip() ):
|
|
|
29 return True
|
|
|
30 else:
|
|
|
31 raise Exception, 'First line in first read does not appear to be a valid FASTQ read in either base-space or color-space'
|
|
|
32 raise Exception, 'There is no non-comment and non-blank line in your FASTQ file'
|
|
|
33
|
|
|
34 def __main__():
|
|
|
35 #Parse Command Line
|
|
|
36 parser = optparse.OptionParser()
|
|
|
37 parser.add_option( '-t', '--threads', dest='threads', help='The number of threads to use' )
|
|
|
38 parser.add_option( '-c', '--color-space', dest='color_space', action='store_true', help='If the input files are SOLiD format' )
|
|
|
39 parser.add_option( '-r', '--ref', dest='ref', help='The reference genome to use or index' )
|
|
|
40 parser.add_option( '-f', '--input1', dest='fastq', help='The (forward) fastq file to use for the mapping' )
|
|
|
41 parser.add_option( '-F', '--input2', dest='rfastq', help='The reverse fastq file to use for mapping if paired-end data' )
|
|
|
42 parser.add_option( '-u', '--output', dest='output', help='The file to save the output (SAM format)' )
|
|
|
43 parser.add_option( '-g', '--genAlignType', dest='genAlignType', help='The type of pairing (single or paired)' )
|
|
|
44 parser.add_option( '-p', '--params', dest='params', help='Parameter setting to use (pre_set or full)' )
|
|
|
45 parser.add_option( '-s', '--fileSource', dest='fileSource', help='Whether to use a previously indexed reference sequence or one form history (indexed or history)' )
|
|
|
46 parser.add_option( '-n', '--maxEditDist', dest='maxEditDist', help='Maximum edit distance if integer' )
|
|
|
47 parser.add_option( '-m', '--fracMissingAligns', dest='fracMissingAligns', help='Fraction of missing alignments given 2% uniform base error rate if fraction' )
|
|
|
48 parser.add_option( '-o', '--maxGapOpens', dest='maxGapOpens', help='Maximum number of gap opens' )
|
|
|
49 parser.add_option( '-e', '--maxGapExtens', dest='maxGapExtens', help='Maximum number of gap extensions' )
|
|
|
50 parser.add_option( '-d', '--disallowLongDel', dest='disallowLongDel', help='Disallow a long deletion within specified bps' )
|
|
|
51 parser.add_option( '-i', '--disallowIndel', dest='disallowIndel', help='Disallow indel within specified bps' )
|
|
|
52 parser.add_option( '-l', '--seed', dest='seed', help='Take the first specified subsequences' )
|
|
|
53 parser.add_option( '-k', '--maxEditDistSeed', dest='maxEditDistSeed', help='Maximum edit distance to the seed' )
|
|
|
54 parser.add_option( '-M', '--mismatchPenalty', dest='mismatchPenalty', help='Mismatch penalty' )
|
|
|
55 parser.add_option( '-O', '--gapOpenPenalty', dest='gapOpenPenalty', help='Gap open penalty' )
|
|
|
56 parser.add_option( '-E', '--gapExtensPenalty', dest='gapExtensPenalty', help='Gap extension penalty' )
|
|
|
57 parser.add_option( '-R', '--suboptAlign', dest='suboptAlign', default=None, help='Proceed with suboptimal alignments even if the top hit is a repeat' )
|
|
|
58 parser.add_option( '-N', '--noIterSearch', dest='noIterSearch', help='Disable iterative search' )
|
|
|
59 parser.add_option( '-T', '--outputTopN', dest='outputTopN', help='Maximum number of alignments to output in the XA tag for reads paired properly' )
|
|
|
60 parser.add_option( '', '--outputTopNDisc', dest='outputTopNDisc', help='Maximum number of alignments to output in the XA tag for disconcordant read pairs (excluding singletons)' )
|
|
|
61 parser.add_option( '-S', '--maxInsertSize', dest='maxInsertSize', help='Maximum insert size for a read pair to be considered mapped good' )
|
|
|
62 parser.add_option( '-P', '--maxOccurPairing', dest='maxOccurPairing', help='Maximum occurrences of a read for pairings' )
|
|
|
63 parser.add_option( '', '--rgid', dest='rgid', help='Read group identifier' )
|
|
|
64 parser.add_option( '', '--rgcn', dest='rgcn', help='Sequencing center that produced the read' )
|
|
|
65 parser.add_option( '', '--rgds', dest='rgds', help='Description' )
|
|
|
66 parser.add_option( '', '--rgdt', dest='rgdt', help='Date that run was produced (ISO8601 format date or date/time, like YYYY-MM-DD)' )
|
|
|
67 parser.add_option( '', '--rgfo', dest='rgfo', help='Flow order' )
|
|
|
68 parser.add_option( '', '--rgks', dest='rgks', help='The array of nucleotide bases that correspond to the key sequence of each read' )
|
|
|
69 parser.add_option( '', '--rglb', dest='rglb', help='Library name' )
|
|
|
70 parser.add_option( '', '--rgpg', dest='rgpg', help='Programs used for processing the read group' )
|
|
|
71 parser.add_option( '', '--rgpi', dest='rgpi', help='Predicted median insert size' )
|
|
|
72 parser.add_option( '', '--rgpl', dest='rgpl', choices=[ 'CAPILLARY', 'LS454', 'ILLUMINA', 'SOLID', 'HELICOS', 'IONTORRENT' and 'PACBIO' ], help='Platform/technology used to produce the reads' )
|
|
|
73 parser.add_option( '', '--rgpu', dest='rgpu', help='Platform unit (e.g. flowcell-barcode.lane for Illumina or slide for SOLiD)' )
|
|
|
74 parser.add_option( '', '--rgsm', dest='rgsm', help='Sample' )
|
|
|
75 parser.add_option( '-D', '--dbkey', dest='dbkey', help='Dbkey for reference genome' )
|
|
|
76 parser.add_option( '-X', '--do_not_build_index', dest='do_not_build_index', action='store_true', help="Don't build index" )
|
|
|
77 parser.add_option( '-H', '--suppressHeader', dest='suppressHeader', help='Suppress header' )
|
|
|
78 parser.add_option( '-I', '--illumina1.3', dest='illumina13qual', help='Input FASTQ files have Illuina 1.3 quality scores' )
|
|
|
79 (options, args) = parser.parse_args()
|
|
|
80
|
|
|
81 # output version # of tool
|
|
|
82 try:
|
|
|
83 tmp = tempfile.NamedTemporaryFile().name
|
|
|
84 tmp_stdout = open( tmp, 'wb' )
|
|
|
85 proc = subprocess.Popen( args='bwa 2>&1', shell=True, stdout=tmp_stdout )
|
|
|
86 tmp_stdout.close()
|
|
|
87 returncode = proc.wait()
|
|
|
88 stdout = None
|
|
|
89 for line in open( tmp_stdout.name, 'rb' ):
|
|
|
90 if line.lower().find( 'version' ) >= 0:
|
|
|
91 stdout = line.strip()
|
|
|
92 break
|
|
|
93 if stdout:
|
|
|
94 sys.stdout.write( 'BWA %s\n' % stdout )
|
|
|
95 else:
|
|
|
96 raise Exception
|
|
|
97 except:
|
|
|
98 sys.stdout.write( 'Could not determine BWA version\n' )
|
|
|
99
|
|
|
100 # check for color space fastq that's not double-encoded and exit if appropriate
|
|
|
101 if options.color_space:
|
|
|
102 if not check_is_double_encoded( options.fastq ):
|
|
|
103 stop_err( 'Your file must be double-encoded (it must be converted from "numbers" to "bases"). See the help section for details' )
|
|
|
104 if options.genAlignType == 'paired':
|
|
|
105 if not check_is_double_encoded( options.rfastq ):
|
|
|
106 stop_err( 'Your reverse reads file must also be double-encoded (it must be converted from "numbers" to "bases"). See the help section for details' )
|
|
|
107
|
|
|
108 fastq = options.fastq
|
|
|
109 if options.rfastq:
|
|
|
110 rfastq = options.rfastq
|
|
|
111
|
|
|
112 # set color space variable
|
|
|
113 if options.color_space:
|
|
|
114 color_space = '-c'
|
|
|
115 else:
|
|
|
116 color_space = ''
|
|
|
117
|
|
|
118 # make temp directory for placement of indices
|
|
|
119 tmp_index_dir = tempfile.mkdtemp()
|
|
|
120 tmp_dir = tempfile.mkdtemp()
|
|
|
121 # index if necessary
|
|
|
122 if options.fileSource == 'history' and not options.do_not_build_index:
|
|
|
123 ref_file = tempfile.NamedTemporaryFile( dir=tmp_index_dir )
|
|
|
124 ref_file_name = ref_file.name
|
|
|
125 ref_file.close()
|
|
|
126 os.symlink( options.ref, ref_file_name )
|
|
|
127 # determine which indexing algorithm to use, based on size
|
|
|
128 try:
|
|
|
129 size = os.stat( options.ref ).st_size
|
|
|
130 if size <= 2**30:
|
|
|
131 indexingAlg = 'is'
|
|
|
132 else:
|
|
|
133 indexingAlg = 'bwtsw'
|
|
|
134 except:
|
|
|
135 indexingAlg = 'is'
|
|
|
136 indexing_cmds = '%s -a %s' % ( color_space, indexingAlg )
|
|
|
137 cmd1 = 'bwa index %s %s' % ( indexing_cmds, ref_file_name )
|
|
|
138 try:
|
|
|
139 tmp = tempfile.NamedTemporaryFile( dir=tmp_index_dir ).name
|
|
|
140 tmp_stderr = open( tmp, 'wb' )
|
|
|
141 proc = subprocess.Popen( args=cmd1, shell=True, cwd=tmp_index_dir, stderr=tmp_stderr.fileno() )
|
|
|
142 returncode = proc.wait()
|
|
|
143 tmp_stderr.close()
|
|
|
144 # get stderr, allowing for case where it's very large
|
|
|
145 tmp_stderr = open( tmp, 'rb' )
|
|
|
146 stderr = ''
|
|
|
147 buffsize = 1048576
|
|
|
148 try:
|
|
|
149 while True:
|
|
|
150 stderr += tmp_stderr.read( buffsize )
|
|
|
151 if not stderr or len( stderr ) % buffsize != 0:
|
|
|
152 break
|
|
|
153 except OverflowError:
|
|
|
154 pass
|
|
|
155 tmp_stderr.close()
|
|
|
156 if returncode != 0:
|
|
|
157 raise Exception, stderr
|
|
|
158 except Exception, e:
|
|
|
159 # clean up temp dirs
|
|
|
160 if os.path.exists( tmp_index_dir ):
|
|
|
161 shutil.rmtree( tmp_index_dir )
|
|
|
162 if os.path.exists( tmp_dir ):
|
|
|
163 shutil.rmtree( tmp_dir )
|
|
|
164 stop_err( 'Error indexing reference sequence. ' + str( e ) )
|
|
|
165 else:
|
|
|
166 ref_file_name = options.ref
|
|
|
167 if options.illumina13qual:
|
|
|
168 illumina_quals = "-I"
|
|
|
169 else:
|
|
|
170 illumina_quals = ""
|
|
|
171
|
|
|
172 # set up aligning and generate aligning command options
|
|
|
173 if options.params == 'pre_set':
|
|
|
174 aligning_cmds = '-t %s %s %s' % ( options.threads, color_space, illumina_quals )
|
|
|
175 gen_alignment_cmds = ''
|
|
|
176 else:
|
|
|
177 if options.maxEditDist != '0':
|
|
|
178 editDist = options.maxEditDist
|
|
|
179 else:
|
|
|
180 editDist = options.fracMissingAligns
|
|
|
181 if options.seed != '-1':
|
|
|
182 seed = '-l %s' % options.seed
|
|
|
183 else:
|
|
|
184 seed = ''
|
|
|
185 if options.suboptAlign:
|
|
|
186 suboptAlign = '-R "%s"' % ( options.suboptAlign )
|
|
|
187 else:
|
|
|
188 suboptAlign = ''
|
|
|
189 if options.noIterSearch == 'true':
|
|
|
190 noIterSearch = '-N'
|
|
|
191 else:
|
|
|
192 noIterSearch = ''
|
|
|
193 aligning_cmds = '-n %s -o %s -e %s -d %s -i %s %s -k %s -t %s -M %s -O %s -E %s %s %s %s %s' % \
|
|
|
194 ( editDist, options.maxGapOpens, options.maxGapExtens, options.disallowLongDel,
|
|
|
195 options.disallowIndel, seed, options.maxEditDistSeed, options.threads,
|
|
|
196 options.mismatchPenalty, options.gapOpenPenalty, options.gapExtensPenalty,
|
|
|
197 suboptAlign, noIterSearch, color_space, illumina_quals )
|
|
|
198 if options.genAlignType == 'paired':
|
|
|
199 gen_alignment_cmds = '-a %s -o %s' % ( options.maxInsertSize, options.maxOccurPairing )
|
|
|
200 if options.outputTopNDisc:
|
|
|
201 gen_alignment_cmds += ' -N %s' % options.outputTopNDisc
|
|
|
202 else:
|
|
|
203 gen_alignment_cmds = ''
|
|
|
204 if options.metadata_tsv:
|
|
|
205 f = open(options.metadata_tsv, 'r')
|
|
|
206 cols = f.readline().split('\t')
|
|
|
207 readGroup = '@RG\tID:1\tDS:%s\tPU:%s\tLB:%s\tCN:BMGC\tDT:%s\tFR:%s\tRR:%s\tPL:ILLUMINA' % ( cols[4], cols[16], cols[14], cols[1], cols[19], cols[20] )
|
|
|
208 gen_alignment_cmds += ' -r "%s"' % readGroup
|
|
|
209 print readGroup
|
|
|
210 if options.rgid:
|
|
|
211 if not options.rglb or not options.rgpl or not options.rgsm:
|
|
|
212 stop_err( 'If you want to specify read groups, you must include the ID, LB, PL, and SM tags.' )
|
|
|
213 readGroup = '@RG\tID:%s\tLB:%s\tPL:%s\tSM:%s' % ( options.rgid, options.rglb, options.rgpl, options.rgsm )
|
|
|
214 if options.rgcn:
|
|
|
215 readGroup += '\tCN:%s' % options.rgcn
|
|
|
216 if options.rgds:
|
|
|
217 readGroup += '\tDS:%s' % options.rgds
|
|
|
218 if options.rgdt:
|
|
|
219 readGroup += '\tDT:%s' % options.rgdt
|
|
|
220 if options.rgfo:
|
|
|
221 readGroup += '\tFO:%s' % options.rgfo
|
|
|
222 if options.rgks:
|
|
|
223 readGroup += '\tKS:%s' % options.rgks
|
|
|
224 if options.rgpg:
|
|
|
225 readGroup += '\tPG:%s' % options.rgpg
|
|
|
226 if options.rgpi:
|
|
|
227 readGroup += '\tPI:%s' % options.rgpi
|
|
|
228 if options.rgpu:
|
|
|
229 readGroup += '\tPU:%s' % options.rgpu
|
|
|
230 gen_alignment_cmds += ' -r "%s"' % readGroup
|
|
|
231 if options.outputTopN:
|
|
|
232 gen_alignment_cmds += ' -n %s' % options.outputTopN
|
|
|
233 # set up output files
|
|
|
234 tmp_align_out = tempfile.NamedTemporaryFile( dir=tmp_dir )
|
|
|
235 tmp_align_out_name = tmp_align_out.name
|
|
|
236 tmp_align_out.close()
|
|
|
237 tmp_align_out2 = tempfile.NamedTemporaryFile( dir=tmp_dir )
|
|
|
238 tmp_align_out2_name = tmp_align_out2.name
|
|
|
239 tmp_align_out2.close()
|
|
|
240 # prepare actual aligning and generate aligning commands
|
|
|
241 cmd2 = 'bwa aln %s %s %s > %s' % ( aligning_cmds, ref_file_name, fastq, tmp_align_out_name )
|
|
|
242 cmd2b = ''
|
|
|
243 if options.genAlignType == 'paired':
|
|
|
244 cmd2b = 'bwa aln %s %s %s > %s' % ( aligning_cmds, ref_file_name, rfastq, tmp_align_out2_name )
|
|
|
245 cmd3 = 'bwa sampe %s %s %s %s %s %s >> %s' % ( gen_alignment_cmds, ref_file_name, tmp_align_out_name, tmp_align_out2_name, fastq, rfastq, options.output )
|
|
|
246 else:
|
|
|
247 cmd3 = 'bwa samse %s %s %s %s >> %s' % ( gen_alignment_cmds, ref_file_name, tmp_align_out_name, fastq, options.output )
|
|
|
248 # perform alignments
|
|
|
249 buffsize = 1048576
|
|
|
250 try:
|
|
|
251 # need to nest try-except in try-finally to handle 2.4
|
|
|
252 try:
|
|
|
253 # align
|
|
|
254 try:
|
|
|
255 tmp = tempfile.NamedTemporaryFile( dir=tmp_dir ).name
|
|
|
256 tmp_stderr = open( tmp, 'wb' )
|
|
|
257 proc = subprocess.Popen( args=cmd2, shell=True, cwd=tmp_dir, stderr=tmp_stderr.fileno() )
|
|
|
258 returncode = proc.wait()
|
|
|
259 tmp_stderr.close()
|
|
|
260 # get stderr, allowing for case where it's very large
|
|
|
261 tmp_stderr = open( tmp, 'rb' )
|
|
|
262 stderr = ''
|
|
|
263 try:
|
|
|
264 while True:
|
|
|
265 stderr += tmp_stderr.read( buffsize )
|
|
|
266 if not stderr or len( stderr ) % buffsize != 0:
|
|
|
267 break
|
|
|
268 except OverflowError:
|
|
|
269 pass
|
|
|
270 tmp_stderr.close()
|
|
|
271 if returncode != 0:
|
|
|
272 raise Exception, stderr
|
|
|
273 except Exception, e:
|
|
|
274 raise Exception, 'Error aligning sequence. ' + str( e )
|
|
|
275 # and again if paired data
|
|
|
276 try:
|
|
|
277 if cmd2b:
|
|
|
278 tmp = tempfile.NamedTemporaryFile( dir=tmp_dir ).name
|
|
|
279 tmp_stderr = open( tmp, 'wb' )
|
|
|
280 proc = subprocess.Popen( args=cmd2b, shell=True, cwd=tmp_dir, stderr=tmp_stderr.fileno() )
|
|
|
281 returncode = proc.wait()
|
|
|
282 tmp_stderr.close()
|
|
|
283 # get stderr, allowing for case where it's very large
|
|
|
284 tmp_stderr = open( tmp, 'rb' )
|
|
|
285 stderr = ''
|
|
|
286 try:
|
|
|
287 while True:
|
|
|
288 stderr += tmp_stderr.read( buffsize )
|
|
|
289 if not stderr or len( stderr ) % buffsize != 0:
|
|
|
290 break
|
|
|
291 except OverflowError:
|
|
|
292 pass
|
|
|
293 tmp_stderr.close()
|
|
|
294 if returncode != 0:
|
|
|
295 raise Exception, stderr
|
|
|
296 except Exception, e:
|
|
|
297 raise Exception, 'Error aligning second sequence. ' + str( e )
|
|
|
298 # generate align
|
|
|
299 try:
|
|
|
300 tmp = tempfile.NamedTemporaryFile( dir=tmp_dir ).name
|
|
|
301 tmp_stderr = open( tmp, 'wb' )
|
|
|
302 proc = subprocess.Popen( args=cmd3, shell=True, cwd=tmp_dir, stderr=tmp_stderr.fileno() )
|
|
|
303 returncode = proc.wait()
|
|
|
304 tmp_stderr.close()
|
|
|
305 # get stderr, allowing for case where it's very large
|
|
|
306 tmp_stderr = open( tmp, 'rb' )
|
|
|
307 stderr = ''
|
|
|
308 try:
|
|
|
309 while True:
|
|
|
310 stderr += tmp_stderr.read( buffsize )
|
|
|
311 if not stderr or len( stderr ) % buffsize != 0:
|
|
|
312 break
|
|
|
313 except OverflowError:
|
|
|
314 pass
|
|
|
315 tmp_stderr.close()
|
|
|
316 if returncode != 0:
|
|
|
317 raise Exception, stderr
|
|
|
318 except Exception, e:
|
|
|
319 raise Exception, 'Error generating alignments. ' + str( e )
|
|
|
320 # remove header if necessary
|
|
|
321 if options.suppressHeader == 'true':
|
|
|
322 tmp_out = tempfile.NamedTemporaryFile( dir=tmp_dir)
|
|
|
323 tmp_out_name = tmp_out.name
|
|
|
324 tmp_out.close()
|
|
|
325 try:
|
|
|
326 shutil.move( options.output, tmp_out_name )
|
|
|
327 except Exception, e:
|
|
|
328 raise Exception, 'Error moving output file before removing headers. ' + str( e )
|
|
|
329 fout = file( options.output, 'w' )
|
|
|
330 for line in file( tmp_out.name, 'r' ):
|
|
|
331 if not ( line.startswith( '@HD' ) or line.startswith( '@SQ' ) or line.startswith( '@RG' ) or line.startswith( '@PG' ) or line.startswith( '@CO' ) ):
|
|
|
332 fout.write( line )
|
|
|
333 fout.close()
|
|
|
334 # check that there are results in the output file
|
|
|
335 if os.path.getsize( options.output ) > 0:
|
|
|
336 sys.stdout.write( 'BWA run on %s-end data' % options.genAlignType )
|
|
|
337 else:
|
|
|
338 raise Exception, 'The output file is empty. You may simply have no matches, or there may be an error with your input file or settings.'
|
|
|
339 except Exception, e:
|
|
|
340 stop_err( 'The alignment failed.\n' + str( e ) )
|
|
|
341 finally:
|
|
|
342 # clean up temp dir
|
|
|
343 if os.path.exists( tmp_index_dir ):
|
|
|
344 shutil.rmtree( tmp_index_dir )
|
|
|
345 if os.path.exists( tmp_dir ):
|
|
|
346 shutil.rmtree( tmp_dir )
|
|
|
347
|
|
|
348 if __name__=="__main__": __main__()
|
