Mercurial > repos > curtish > bwa_wrappers_demo

comparison bwa_wrapper.py @ 0:cdf2720a3d62 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Prepare to Fork off main branch
author curtish
date Mon, 26 Aug 2013 16:49:03 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:cdf2720a3d62
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.rgid:
205 if not options.rglb or not options.rgpl or not options.rgsm:
206 stop_err( 'If you want to specify read groups, you must include the ID, LB, PL, and SM tags.' )
207 readGroup = '@RG\tID:%s\tLB:%s\tPL:%s\tSM:%s' % ( options.rgid, options.rglb, options.rgpl, options.rgsm )
208 if options.rgcn:
209 readGroup += '\tCN:%s' % options.rgcn
210 if options.rgds:
211 readGroup += '\tDS:%s' % options.rgds
212 if options.rgdt:
213 readGroup += '\tDT:%s' % options.rgdt
214 if options.rgfo:
215 readGroup += '\tFO:%s' % options.rgfo
216 if options.rgks:
217 readGroup += '\tKS:%s' % options.rgks
218 if options.rgpg:
219 readGroup += '\tPG:%s' % options.rgpg
220 if options.rgpi:
221 readGroup += '\tPI:%s' % options.rgpi
222 if options.rgpu:
223 readGroup += '\tPU:%s' % options.rgpu
224 gen_alignment_cmds += ' -r "%s"' % readGroup
225 if options.outputTopN:
226 gen_alignment_cmds += ' -n %s' % options.outputTopN
227 # set up output files
228 tmp_align_out = tempfile.NamedTemporaryFile( dir=tmp_dir )
229 tmp_align_out_name = tmp_align_out.name
230 tmp_align_out.close()
231 tmp_align_out2 = tempfile.NamedTemporaryFile( dir=tmp_dir )
232 tmp_align_out2_name = tmp_align_out2.name
233 tmp_align_out2.close()
234 # prepare actual aligning and generate aligning commands
235 cmd2 = 'bwa aln %s %s %s > %s' % ( aligning_cmds, ref_file_name, fastq, tmp_align_out_name )
236 cmd2b = ''
237 if options.genAlignType == 'paired':
238 cmd2b = 'bwa aln %s %s %s > %s' % ( aligning_cmds, ref_file_name, rfastq, tmp_align_out2_name )
239 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 )
240 else:
241 cmd3 = 'bwa samse %s %s %s %s >> %s' % ( gen_alignment_cmds, ref_file_name, tmp_align_out_name, fastq, options.output )
242 # perform alignments
243 buffsize = 1048576
244 try:
245 # need to nest try-except in try-finally to handle 2.4
246 try:
247 # align
248 try:
249 tmp = tempfile.NamedTemporaryFile( dir=tmp_dir ).name
250 tmp_stderr = open( tmp, 'wb' )
251 proc = subprocess.Popen( args=cmd2, shell=True, cwd=tmp_dir, stderr=tmp_stderr.fileno() )
252 returncode = proc.wait()
253 tmp_stderr.close()
254 # get stderr, allowing for case where it's very large
255 tmp_stderr = open( tmp, 'rb' )
256 stderr = ''
257 try:
258 while True:
259 stderr += tmp_stderr.read( buffsize )
260 if not stderr or len( stderr ) % buffsize != 0:
261 break
262 except OverflowError:
263 pass
264 tmp_stderr.close()
265 if returncode != 0:
266 raise Exception, stderr
267 except Exception, e:
268 raise Exception, 'Error aligning sequence. ' + str( e )
269 # and again if paired data
270 try:
271 if cmd2b:
272 tmp = tempfile.NamedTemporaryFile( dir=tmp_dir ).name
273 tmp_stderr = open( tmp, 'wb' )
274 proc = subprocess.Popen( args=cmd2b, shell=True, cwd=tmp_dir, stderr=tmp_stderr.fileno() )
275 returncode = proc.wait()
276 tmp_stderr.close()
277 # get stderr, allowing for case where it's very large
278 tmp_stderr = open( tmp, 'rb' )
279 stderr = ''
280 try:
281 while True:
282 stderr += tmp_stderr.read( buffsize )
283 if not stderr or len( stderr ) % buffsize != 0:
284 break
285 except OverflowError:
286 pass
287 tmp_stderr.close()
288 if returncode != 0:
289 raise Exception, stderr
290 except Exception, e:
291 raise Exception, 'Error aligning second sequence. ' + str( e )
292 # generate align
293 try:
294 tmp = tempfile.NamedTemporaryFile( dir=tmp_dir ).name
295 tmp_stderr = open( tmp, 'wb' )
296 proc = subprocess.Popen( args=cmd3, shell=True, cwd=tmp_dir, stderr=tmp_stderr.fileno() )
297 returncode = proc.wait()
298 tmp_stderr.close()
299 # get stderr, allowing for case where it's very large
300 tmp_stderr = open( tmp, 'rb' )
301 stderr = ''
302 try:
303 while True:
304 stderr += tmp_stderr.read( buffsize )
305 if not stderr or len( stderr ) % buffsize != 0:
306 break
307 except OverflowError:
308 pass
309 tmp_stderr.close()
310 if returncode != 0:
311 raise Exception, stderr
312 except Exception, e:
313 raise Exception, 'Error generating alignments. ' + str( e )
314 # remove header if necessary
315 if options.suppressHeader == 'true':
316 tmp_out = tempfile.NamedTemporaryFile( dir=tmp_dir)
317 tmp_out_name = tmp_out.name
318 tmp_out.close()
319 try:
320 shutil.move( options.output, tmp_out_name )
321 except Exception, e:
322 raise Exception, 'Error moving output file before removing headers. ' + str( e )
323 fout = file( options.output, 'w' )
324 for line in file( tmp_out.name, 'r' ):
325 if not ( line.startswith( '@HD' ) or line.startswith( '@SQ' ) or line.startswith( '@RG' ) or line.startswith( '@PG' ) or line.startswith( '@CO' ) ):
326 fout.write( line )
327 fout.close()
328 # check that there are results in the output file
329 if os.path.getsize( options.output ) > 0:
330 sys.stdout.write( 'BWA run on %s-end data' % options.genAlignType )
331 else:
332 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.'
333 except Exception, e:
334 stop_err( 'The alignment failed.\n' + str( e ) )
335 finally:
336 # clean up temp dir
337 if os.path.exists( tmp_index_dir ):
338 shutil.rmtree( tmp_index_dir )
339 if os.path.exists( tmp_dir ):
340 shutil.rmtree( tmp_dir )
341
342 if __name__=="__main__": __main__()