Mercurial > repos > crs4 > bwa_mem

comparison bwa_mem.py @ 0:6820983ba5d5 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author crs4
date Tue, 18 Mar 2014 07:49:22 -0400
parents
children ebb02ba5987c
comparison
equal deleted inserted replaced
-1:000000000000 0:6820983ba5d5
1 # -*- coding: utf-8 -*-
2 #!/usr/bin/env python
3 ## yufei.luo@gustave.roussy 22/07/2013
4 ## Copyright © 2014 CRS4 Srl. http://www.crs4.it/
5 ## Modified by:
6 ## Nicola Soranzo <nicola.soranzo@crs4.it>
7
8 """
9 Runs BWA on single-end or paired-end data.
10 Produces a SAM file containing the mappings.
11 Works with BWA version >= 0.7.5.
12 NOTICE: In this wrapper, we only use 'mem' for mapping step.
13
14 usage: bwa_mem.py [options]
15
16 See below for options
17 """
18
19 import optparse, os, shutil, subprocess, sys, tempfile
20
21 def stop_err( msg ):
22 sys.stderr.write( '%s\n' % msg )
23 sys.exit()
24
25 def check_is_double_encoded( fastq ):
26 # check that first read is bases, not one base followed by numbers
27 bases = [ 'A', 'C', 'G', 'T', 'a', 'c', 'g', 't', 'N' ]
28 nums = [ '0', '1', '2', '3' ]
29 for line in file( fastq, 'rb'):
30 if not line.strip() or line.startswith( '@' ):
31 continue
32 if len( [ b for b in line.strip() if b in nums ] ) > 0:
33 return False
34 elif line.strip()[0] in bases and len( [ b for b in line.strip() if b in bases ] ) == len( line.strip() ):
35 return True
36 else:
37 raise Exception, 'First line in first read does not appear to be a valid FASTQ read in either base-space or color-space'
38 raise Exception, 'There is no non-comment and non-blank line in your FASTQ file'
39
40 def __main__():
41 descr = "bwa_mem.py: Map (long length) reads against a reference genome with BWA-MEM."
42 parser = optparse.OptionParser(description=descr)
43 parser.add_option( '-t', '--threads', default=1, help='The number of threads to use [1]' )
44 parser.add_option( '--ref', help='The reference genome to use or index' )
45 parser.add_option( '-f', '--fastq', help='The (forward) fastq file to use for the mapping' )
46 parser.add_option( '-F', '--rfastq', help='The reverse fastq file to use for mapping if paired-end data' )
47 parser.add_option( '-u', '--output', help='The file to save the output (SAM format)' )
48 parser.add_option( '-g', '--genAlignType', help='The type of pairing (single or paired)' )
49 parser.add_option( '--params', help='Parameter setting to use (pre_set or full)' )
50 parser.add_option( '-s', '--fileSource', help='Whether to use a previously indexed reference sequence or one form history (indexed or history)' )
51 parser.add_option( '-D', '--dbkey', help='Dbkey for reference genome' )
52
53 parser.add_option( '-k', '--minEditDistSeed', default=19, type=int, help='Minimum edit distance to the seed [19]' )
54 parser.add_option( '-w', '--bandWidth', default=100, type=int, help='Band width for banded alignment [100]' )
55 parser.add_option( '-d', '--offDiagonal', default=100, type=int, help='off-diagonal X-dropoff [100]' )
56 parser.add_option( '-r', '--internalSeeds', default=1.5, type=float, help='look for internal seeds inside a seed longer than {-k} * FLOAT [1.5]' )
57 parser.add_option( '-c', '--seedsOccurrence', default=10000, type=int, help='skip seeds with more than INT occurrences [10000]' )
58 parser.add_option( '-S', '--mateRescue', default=False, help='skip mate rescue' )
59 parser.add_option( '-P', '--skipPairing', default=False, help='skpe pairing, mate rescue performed unless -S also in use' )
60 parser.add_option( '-A', '--seqMatch', default=1, type=int, help='score of a sequence match' )
61 parser.add_option( '-B', '--mismatch', default=4, type=int, help='penalty for a mismatch' )
62 parser.add_option( '-O', '--gapOpen', default=6, type=int, help='gap open penalty' )
63 parser.add_option( '-E', '--gapExtension', default=None, help='gap extension penalty; a gap of size k cost {-O} + {-E}*k [1]' )
64 parser.add_option( '-L', '--clipping', default=5, type=int, help='penalty for clipping [5]' )
65 parser.add_option( '-U', '--unpairedReadpair', default=17, type=int, help='penalty for an unpaired read pair [17]' )
66 parser.add_option( '-p', '--interPairEnd', default=False, help='first query file consists of interleaved paired-end sequences' )
67 parser.add_option( '--rgid', help='Read group identifier' )
68 parser.add_option( '--rgsm', help='Sample' )
69 parser.add_option( '--rgpl', choices=[ 'CAPILLARY', 'LS454', 'ILLUMINA', 'SOLID', 'HELICOS', 'IONTORRENT', 'PACBIO' ], help='Platform/technology used to produce the reads' )
70 parser.add_option( '--rglb', help='Library name' )
71 parser.add_option( '--rgpu', help='Platform unit (e.g. flowcell-barcode.lane for Illumina or slide for SOLiD)' )
72 parser.add_option( '--rgcn', help='Sequencing center that produced the read' )
73 parser.add_option( '--rgds', help='Description' )
74 parser.add_option( '--rgdt', help='Date that run was produced (ISO8601 format date or date/time, like YYYY-MM-DD)' )
75 parser.add_option( '--rgfo', help='Flow order' )
76 parser.add_option( '--rgks', help='The array of nucleotide bases that correspond to the key sequence of each read' )
77 parser.add_option( '--rgpg', help='Programs used for processing the read group' )
78 parser.add_option( '--rgpi', help='Predicted median insert size' )
79 parser.add_option( '-T', '--minScore', default=30, type=int, help='minimum score to output [30]' )
80 parser.add_option( '-M', '--mark', default=False, help='mark shorter split hits as secondary (for Picard/GATK compatibility)' )
81 parser.add_option( '-H', '--suppressHeader', dest='suppressHeader', help='Suppress header' )
82 (options, args) = parser.parse_args()
83 if len(args) > 0:
84 parser.error('Wrong number of arguments')
85
86 # output version # of tool
87 try:
88 tmp = tempfile.NamedTemporaryFile().name
89 tmp_stdout = open( tmp, 'wb' )
90 proc = subprocess.Popen( args='bwa 2>&1', shell=True, stdout=tmp_stdout )
91 tmp_stdout.close()
92 returncode = proc.wait()
93 stdout = None
94 for line in open( tmp_stdout.name, 'rb' ):
95 if line.lower().find( 'version' ) >= 0:
96 stdout = line.strip()
97 break
98 if stdout:
99 sys.stdout.write( 'BWA %s\n' % stdout )
100 else:
101 raise Exception
102 except:
103 sys.stdout.write( 'Could not determine BWA version\n' )
104
105 fastq = options.fastq
106 if options.rfastq:
107 rfastq = options.rfastq
108
109 # make temp directory for placement of indices
110 tmp_index_dir = tempfile.mkdtemp()
111 tmp_dir = tempfile.mkdtemp()
112 # index if necessary
113 if options.fileSource == 'history':
114 ref_file = tempfile.NamedTemporaryFile( dir=tmp_index_dir )
115 ref_file_name = ref_file.name
116 ref_file.close()
117 os.symlink( options.ref, ref_file_name )
118 # determine which indexing algorithm to use, based on size
119 try:
120 size = os.stat( options.ref ).st_size
121 if size <= 2**30:
122 indexingAlg = 'is'
123 else:
124 indexingAlg = 'bwtsw'
125 except:
126 indexingAlg = 'is'
127 indexing_cmds = '-a %s' % indexingAlg
128 cmd1 = 'bwa index %s %s' % ( indexing_cmds, ref_file_name )
129 try:
130 tmp = tempfile.NamedTemporaryFile( dir=tmp_index_dir ).name
131 tmp_stderr = open( tmp, 'wb' )
132 proc = subprocess.Popen( args=cmd1, shell=True, cwd=tmp_index_dir, stderr=tmp_stderr.fileno() )
133 returncode = proc.wait()
134 tmp_stderr.close()
135 # get stderr, allowing for case where it's very large
136 tmp_stderr = open( tmp, 'rb' )
137 stderr = ''
138 buffsize = 1048576
139 try:
140 while True:
141 stderr += tmp_stderr.read( buffsize )
142 if not stderr or len( stderr ) % buffsize != 0:
143 break
144 except OverflowError:
145 pass
146 tmp_stderr.close()
147 if returncode != 0:
148 raise Exception, stderr
149 except Exception, e:
150 # clean up temp dirs
151 if os.path.exists( tmp_index_dir ):
152 shutil.rmtree( tmp_index_dir )
153 if os.path.exists( tmp_dir ):
154 shutil.rmtree( tmp_dir )
155 stop_err( 'Error indexing reference sequence. ' + str( e ) )
156 else:
157 ref_file_name = options.ref
158 # if options.illumina13qual:
159 # illumina_quals = "-I"
160 # else:
161 # illumina_quals = ""
162
163 # set up aligning and generate aligning command args
164 start_cmds = '-t %s ' % options.threads
165 if options.params == 'pre_set':
166 # aligning_cmds = '-t %s %s' % ( options.threads, illumina_quals )
167 #start_cmds = '-t %s ' % options.threads
168 end_cmds = ' '
169 print start_cmds, end_cmds
170
171 else:
172 end_cmds = '-k %s -w %s -d %s -r %s -c %s -A %s -B %s -O %s -L %s -U %s -T %s ' % (options.minEditDistSeed, options.bandWidth, options.offDiagonal, options.internalSeeds, options.seedsOccurrence, options.seqMatch, options.mismatch, options.gapOpen, options.clipping, options.unpairedReadpair, options.minScore)
173 if options.mateRescue:
174 end_cmds += '-S '
175 if options.skipPairing:
176 end_cmds += '-P '
177 else:
178 if options.skipPairing:
179 print "Option Error and will not be considered, you should also choose 'skip mate rescue -S' option! "
180 if options.gapExtension != None:
181 end_cmds += '-E %s ' % options.gapExtension
182
183 if options.rgid:
184 if not options.rglb or not options.rgpl or not options.rgsm or not options.rglb:
185 stop_err( 'If you want to specify read groups, you must include the ID, LB, PL, and SM tags.' )
186 # readGroup = '@RG\tID:%s\tLB:%s\tPL:%s\tSM:%s' % ( options.rgid, options.rglb, options.rgpl, options.rgsm )
187 readGroup = '@RG\tID:%s\tLB:%s\tPL:%s\tSM:%s' % ( options.rgid, options.rglb, options.rgpl, options.rgsm )
188 if options.rgpu:
189 readGroup += '\tPU:%s' % options.rgpu
190 if options.rgcn:
191 readGroup += '\tCN:%s' % options.rgcn
192 if options.rgds:
193 readGroup += '\tDS:%s' % options.rgds
194 if options.rgdt:
195 readGroup += '\tDT:%s' % options.rgdt
196 if options.rgfo:
197 readGroup += '\tFO:%s' % options.rgfo
198 if options.rgks:
199 readGroup += '\tKS:%s' % options.rgks
200 if options.rgpg:
201 readGroup += '\tPG:%s' % options.rgpg
202 if options.rgpi:
203 readGroup += '\tPI:%s' % options.rgpi
204 end_cmds += ' -R "%s" ' % readGroup
205
206 if options.interPairEnd:
207 end_cmds += '-p %s ' % options.interPairEnd
208 if options.mark:
209 end_cmds += '-M '
210
211
212 if options.genAlignType == 'paired':
213 cmd = 'bwa mem %s %s %s %s %s > %s' % ( start_cmds, ref_file_name, fastq, rfastq, end_cmds, options.output )
214 else:
215 cmd = 'bwa mem %s %s %s %s > %s' % ( start_cmds, ref_file_name, fastq, end_cmds, options.output )
216
217 # perform alignments
218 buffsize = 1048576
219 try:
220 # need to nest try-except in try-finally to handle 2.4
221 try:
222 try:
223 tmp = tempfile.NamedTemporaryFile( dir=tmp_dir ).name
224 tmp_stderr = open( tmp, 'wb' )
225 print "The cmd is %s" % cmd
226 proc = subprocess.Popen( args=cmd, shell=True, cwd=tmp_dir, stderr=tmp_stderr.fileno() )
227 returncode = proc.wait()
228 tmp_stderr.close()
229 # get stderr, allowing for case where it's very large
230 tmp_stderr = open( tmp, 'rb' )
231 stderr = ''
232 try:
233 while True:
234 stderr += tmp_stderr.read( buffsize )
235 if not stderr or len( stderr ) % buffsize != 0:
236 break
237 except OverflowError:
238 pass
239 tmp_stderr.close()
240 if returncode != 0:
241 raise Exception, stderr
242 except Exception, e:
243 raise Exception, 'Error generating alignments. ' + str( e )
244 # remove header if necessary
245 if options.suppressHeader == 'true':
246 tmp_out = tempfile.NamedTemporaryFile( dir=tmp_dir)
247 tmp_out_name = tmp_out.name
248 tmp_out.close()
249 try:
250 shutil.move( options.output, tmp_out_name )
251 except Exception, e:
252 raise Exception, 'Error moving output file before removing headers. ' + str( e )
253 fout = file( options.output, 'w' )
254 for line in file( tmp_out.name, 'r' ):
255 if not ( line.startswith( '@HD' ) or line.startswith( '@SQ' ) or line.startswith( '@RG' ) or line.startswith( '@PG' ) or line.startswith( '@CO' ) ):
256 fout.write( line )
257 fout.close()
258 # check that there are results in the output file
259 if os.path.getsize( options.output ) > 0:
260 sys.stdout.write( 'BWA run on %s-end data' % options.genAlignType )
261 else:
262 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.'
263 except Exception, e:
264 stop_err( 'The alignment failed.\n' + str( e ) )
265 finally:
266 # clean up temp dir
267 if os.path.exists( tmp_index_dir ):
268 shutil.rmtree( tmp_index_dir )
269 if os.path.exists( tmp_dir ):
270 shutil.rmtree( tmp_dir )
271
272 if __name__ == "__main__":
273 __main__()