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