Mercurial > repos > yufei-luo > bwa_0_7_5
annotate bwa_0_7_5/bwa_0_7_5.py @ 4:22449c3a0b7f draft default tip
debug when use 'history' input fasta as genome reference
author | yufei-luo |
---|---|
date | Tue, 10 Dec 2013 08:48:23 -0500 |
parents | 8409cff2d740 |
children |
rev | line source |
---|---|
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' ) | |
4
22449c3a0b7f
debug when use 'history' input fasta as genome reference
yufei-luo
parents:
3
diff
changeset
|
69 parser.add_argument( '--rgpl', choices=[ 'CAPILLARY', 'LS454', 'ILLUMINA', 'SOLID', 'HELICOS', 'IONTORRENT', 'PACBIO' ], help='Platform/technology used to produce the reads' ) |
0 | 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 | |
4
22449c3a0b7f
debug when use 'history' input fasta as genome reference
yufei-luo
parents:
3
diff
changeset
|
125 # if args.fileSource == 'history' and not args.do_not_build_index: |
22449c3a0b7f
debug when use 'history' input fasta as genome reference
yufei-luo
parents:
3
diff
changeset
|
126 if args.fileSource == 'history' : |
0 | 127 ref_file = tempfile.NamedTemporaryFile( dir=tmp_index_dir ) |
128 ref_file_name = ref_file.name | |
129 ref_file.close() | |
130 os.symlink( args.ref, ref_file_name ) | |
131 # determine which indexing algorithm to use, based on size | |
132 try: | |
133 size = os.stat( args.ref ).st_size | |
134 if size <= 2**30: | |
135 indexingAlg = 'is' | |
136 else: | |
137 indexingAlg = 'bwtsw' | |
138 except: | |
139 indexingAlg = 'is' | |
140 #indexing_cmds = '%s -a %s' % ( color_space, indexingAlg ) | |
141 indexing_cmds = '-a %s' % indexingAlg | |
142 cmd1 = 'bwa index %s %s' % ( indexing_cmds, ref_file_name ) | |
143 try: | |
144 tmp = tempfile.NamedTemporaryFile( dir=tmp_index_dir ).name | |
145 tmp_stderr = open( tmp, 'wb' ) | |
146 proc = subprocess.Popen( args=cmd1, shell=True, cwd=tmp_index_dir, stderr=tmp_stderr.fileno() ) | |
147 returncode = proc.wait() | |
148 tmp_stderr.close() | |
149 # get stderr, allowing for case where it's very large | |
150 tmp_stderr = open( tmp, 'rb' ) | |
151 stderr = '' | |
152 buffsize = 1048576 | |
153 try: | |
154 while True: | |
155 stderr += tmp_stderr.read( buffsize ) | |
156 if not stderr or len( stderr ) % buffsize != 0: | |
157 break | |
158 except OverflowError: | |
159 pass | |
160 tmp_stderr.close() | |
161 if returncode != 0: | |
162 raise Exception, stderr | |
163 except Exception, e: | |
164 # clean up temp dirs | |
165 if os.path.exists( tmp_index_dir ): | |
166 shutil.rmtree( tmp_index_dir ) | |
167 if os.path.exists( tmp_dir ): | |
168 shutil.rmtree( tmp_dir ) | |
169 stop_err( 'Error indexing reference sequence. ' + str( e ) ) | |
170 else: | |
171 ref_file_name = args.ref | |
172 # if args.illumina13qual: | |
173 # illumina_quals = "-I" | |
174 # else: | |
175 # illumina_quals = "" | |
176 | |
177 # set up aligning and generate aligning command args | |
178 start_cmds = '-t %s ' % args.threads | |
179 if args.params == 'pre_set': | |
180 # aligning_cmds = '-t %s %s %s' % ( args.threads, color_space, illumina_quals ) | |
181 #start_cmds = '-t %s ' % args.threads | |
182 end_cmds = ' ' | |
183 print start_cmds, end_cmds | |
184 | |
185 else: | |
186 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) | |
187 if args.mateRescue: | |
188 end_cmds += '-S ' | |
189 if args.skipPairing: | |
190 end_cmds += '-P ' | |
191 else: | |
192 if args.skipPairing: | |
193 print "Option Error and will not be considered, you should also choose 'skip mate rescue -S' option! " | |
194 if args.gapExtension != None: | |
195 end_cmds += '-E %s ' % args.gapExtension | |
196 | |
197 if args.rgid: | |
198 if not args.rglb or not args.rgpl or not args.rgsm or not args.rglb: | |
199 stop_err( 'If you want to specify read groups, you must include the ID, LB, PL, and SM tags.' ) | |
200 # readGroup = '@RG\tID:%s\tLB:%s\tPL:%s\tSM:%s' % ( args.rgid, args.rglb, args.rgpl, args.rgsm ) | |
201 readGroup = '@RG\tID:%s\tLB:%s\tPL:%s\tSM:%s' % ( args.rgid, args.rglb, args.rgpl, args.rgsm ) | |
202 if args.rgpu: | |
203 readGroup += '\tPU:%s' % args.rgpu | |
204 if args.rgcn: | |
205 readGroup += '\tCN:%s' % args.rgcn | |
206 if args.rgds: | |
207 readGroup += '\tDS:%s' % args.rgds | |
208 if args.rgdt: | |
209 readGroup += '\tDT:%s' % args.rgdt | |
210 if args.rgfo: | |
211 readGroup += '\tFO:%s' % args.rgfo | |
212 if args.rgks: | |
213 readGroup += '\tKS:%s' % args.rgks | |
214 if args.rgpg: | |
215 readGroup += '\tPG:%s' % args.rgpg | |
216 if args.rgpi: | |
217 readGroup += '\tPI:%s' % args.rgpi | |
218 end_cmds += ' -R "%s" ' % readGroup | |
219 | |
220 if args.interPairEnd: | |
221 end_cmds += '-p %s ' % args.interPairEnd | |
222 if args.mark: | |
223 end_cmds += '-M ' | |
224 | |
225 | |
226 if args.genAlignType == 'paired': | |
227 cmd = 'bwa mem %s %s %s %s %s > %s' % ( start_cmds, ref_file_name, fastq, rfastq, end_cmds, args.output ) | |
228 else: | |
3 | 229 cmd = 'bwa mem %s %s %s %s > %s' % ( start_cmds, ref_file_name, fastq, end_cmds, args.output ) |
0 | 230 |
231 # perform alignments | |
232 buffsize = 1048576 | |
233 try: | |
234 # need to nest try-except in try-finally to handle 2.4 | |
235 try: | |
236 try: | |
237 tmp = tempfile.NamedTemporaryFile( dir=tmp_dir ).name | |
238 tmp_stderr = open( tmp, 'wb' ) | |
239 print "The cmd is %s" % cmd | |
240 proc = subprocess.Popen( args=cmd, shell=True, cwd=tmp_dir, stderr=tmp_stderr.fileno() ) | |
241 returncode = proc.wait() | |
242 tmp_stderr.close() | |
243 # get stderr, allowing for case where it's very large | |
244 tmp_stderr = open( tmp, 'rb' ) | |
245 stderr = '' | |
246 try: | |
247 while True: | |
248 stderr += tmp_stderr.read( buffsize ) | |
249 if not stderr or len( stderr ) % buffsize != 0: | |
250 break | |
251 except OverflowError: | |
252 pass | |
253 tmp_stderr.close() | |
254 if returncode != 0: | |
255 raise Exception, stderr | |
256 except Exception, e: | |
257 raise Exception, 'Error generating alignments. ' + str( e ) | |
258 | |
259 # check that there are results in the output file | |
260 if os.path.getsize( args.output ) > 0: | |
261 sys.stdout.write( 'BWA run on %s-end data' % args.genAlignType ) | |
262 else: | |
263 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.' | |
264 except Exception, e: | |
265 stop_err( 'The alignment failed.\n' + str( e ) ) | |
266 finally: | |
267 # clean up temp dir | |
268 if os.path.exists( tmp_index_dir ): | |
269 shutil.rmtree( tmp_index_dir ) | |
270 if os.path.exists( tmp_dir ): | |
271 shutil.rmtree( tmp_dir ) | |
272 | |
273 if __name__=="__main__": __main__() |