Mercurial > repos > nilshomer > tmap_wrapper

comparison tmap_wrapper_0.0.19/tmap_wrapper.py @ 0:7cbfd271d207

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Migrated tool version 0.0.19 from old tool shed archive to new tool shed repository
author nilshomer
date Tue, 07 Jun 2011 17:28:45 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:7cbfd271d207
1 #!/usr/bin/env python
2
3 # TODO
4 # - map1/map2/map3 specific options
5
6 """
7 Runs TMAP on Ion Torrent data.
8 Produces a SAM file containing the mappings.
9 Works with TMAP version 0.0.19 or higher.
10
11 usage: tmap_wrapper.py [options]
12 --threads: The number of threads to use
13 --ref: The reference genome to use or index
14 --fastq: The fastq file to use for the mapping
15 --output: The file to save the output (SAM format)
16 --params: Parameter setting to use (pre_set or full)
17 --fileSource: Whether to use a previously indexed reference sequence or one from history (indexed or history)
18 --matchScore: The match score
19 --mismatchPenalty: The ismatch penalty
20 --gapOpenPenalty: The gap open penalty
21 --gapExtensPenalty: The gap extension penalty
22 --flowPenalty: The flow score penalty
23 --flowOrder: The flow order ([ACGT]{4+} or "sff")
24 --bandWidth: The band width
25 --globalMap: The soft-clipping type (0 - allow on the right and left, 1 - allow on the left, 2 - allow on the right, 3 - do not allow soft-clipping)
26 --duplicateWindow: Remove duplicate alignments from different algorithms within this bp window (-1 to disable)
27 --scoringThreshold: The score threshold divided by the match score
28 --queueSize: The queue size for the reads
29 --outputFilter: The output filter (0 - unique best hits, 1 - random best hit, 2 - all best htis, 3 - all alignments)
30 --rgTag: The flag to specify RG tag(s) in the SAM header
31 --rgTagID: The RG ID tag to add to the SAM header
32 --rgTagCN: The RG CN tag to add to the SAM header
33 --rgTagDS: The RG DS tag to add to the SAM header
34 --rgTagDT: The RG DT tag to add to the SAM header
35 --rgTagLB: The RG LB tag to add to the SAM header
36 --rgTagPI: The RG PI tag to add to the SAM header
37 --rgTagPL: The RG PL tag to add to the SAM header
38 --rgTagPU: The RG PU tag to add to the SAM header
39 --rgTagSM: The RG SM tag to add to the SAM header
40 --filterIndependently: Apply the output filter for each algorithm independently
41
42 --map1: Flag to run map1 in the first stage
43 --map1SeedLength: The k-mer length to seed CALs (-1 to disable)
44 --map1SeedMismatches: The maximum number of mismatches in the seed
45 --map1SecondarySeedLength: The secondary seed length (-1 to disable)
46 --map1NumEdits: The maximum number of edits or false-negative probability assuming the maximum error rate
47 --map1BaseError: The assumed per-base maximum error rate
48 --map1Mismatches: The maximum number of or (read length) fraction of mismatches
49 --map1GapOpens: The maximum number of or (read length) fraction of indel starts
50 --map1GapExtensions: The maximum number of or (read length) fraction of indel extensions
51 --map1MaxCALsDeletion: The maximum number of CALs to extend a deletion
52 --map1EndIndels: Indels are not allowed within this number of bps from the end of the read
53 --map1MaxOptimalCALs: Stop searching when INT optimal CALs have been found
54 --map1MaxNodes: The maximum number of alignment nodes
55
56 --map2: Flag to run map2 in the first stage
57 --map2Coefficient: The coefficient of length-threshold adjustment
58 --map2SeedIntervalSize: The maximum seeding interval size
59 --map2ZBest: Keep the z-best nodes during prefix trie traversal
60 --map2ReverseTrigger: The # seeds to trigger reverse alignment
61
62 --map3: Flag to run map3 in the first stage
63 --map3SeedLength: The k-mer length to seed CALs (-1 tunes to the genome size)
64 --map3SeedMaxHits: The maximum number of hits returned by a seed
65 --map3SeedWindow: The window of bases in which to group seeds
66 --map3HPEnumeration: The single homopolymer error difference for enumeration
67
68 --MAP1: Flag to run MAP1 in the second stage
69 --MAP1SeedLength: The k-mer length to seed CALs (-1 to disable)
70 --MAP1SeedMismatches: The maximum number of mismatches in the seed
71 --MAP1SecondarySeedLength: The secondary seed length (-1 to disable)
72 --MAP1NumEdits: The maximum number of edits or false-negative probability assuming the maximum error rate
73 --MAP1BaseError: The assumed per-base maximum error rate
74 --MAP1Mismatches: The maximum number of or (read length) fraction of mismatches
75 --MAP1GapOpens: The maximum number of or (read length) fraction of indel starts
76 --MAP1GapExtensions: The maximum number of or (read length) fraction of indel extensions
77 --MAP1MaxCALsDeletion: The maximum number of CALs to extend a deletion
78 --MAP1EndIndels: Indels are not allowed within this number of bps from the end of the read
79 --MAP1MaxOptimalCALs: Stop searching when INT optimal CALs have been found
80 --MAP1MaxNodes: The maximum number of alignment nodes
81
82 --MAP2: Flag to run MAP2 in the second stage
83 --MAP2Coefficient: The coefficient of length-threshold adjustment
84 --MAP2SeedIntervalSize: The maximum seeding interval size
85 --MAP2ZBest: Keep the z-best nodes during prefix trie traversal
86 --MAP2ReverseTrigger: The # seeds to trigger reverse alignment
87
88 --MAP3: Flag to run MAP3 in the second stage
89 --MAP3SeedLength: The k-mer length to seed CALs (-1 tunes to the genome size)
90 --MAP3SeedMaxHits: The maximum number of hits returned by a seed
91 --MAP3SeedWindow: The window of bases in which to group seeds
92 --MAP3HPEnumeration: The single homopolymer error difference for enumeration
93
94 --suppressHeader: Suppress header
95 --dbkey: Dbkey for reference genome
96 --do_not_build_index: Flag to specify that provided file is already indexed and to just use 'as is'
97 """
98
99 import optparse, os, shutil, subprocess, sys, tempfile
100
101 def stop_err( msg ):
102 sys.stderr.write( '%s\n' % msg )
103 sys.exit()
104
105 def map1_parse( map1SeedLength, map1SeedMismatches, \
106 map1SecondarySeedLength, map1NumEdits, \
107 map1BaseError, map1Mismatches, \
108 map1GapOpens, map1GapExtensions, map1MaxCALsDeletion, \
109 map1EndIndels, map1MaxOptimalCALs, map1MaxNodes ):
110 return '-l %s -s %s -L %s -p %s -P %s -m %s -o %s -e %s -d %s -i %s -b %s -Q %s' % \
111 (map1SeedLength, map1SeedMismatches, \
112 map1SecondarySeedLength, map1NumEdits, \
113 map1BaseError, map1Mismatches, \
114 map1GapOpens, map1GapExtensions, map1MaxCALsDeletion, \
115 map1EndIndels, map1MaxOptimalCALs, map1MaxNodes )
116
117 def map2_parse( map2Coefficient, map2SeedIntervalSize, map2ZBest, map2ReverseTrigger ):
118 return '-c %s -S %s -b %s -N %s' % ( map2Coefficient, map2SeedIntervalSize, map2ZBest, map2ReverseTrigger )
119
120 def map3_parse( map3SeedLength, map3SeedMaxHits, map3SeedWindow, map3HPEnumeration ):
121 return '-l %s -S %s -b %s -H %s' % ( map3SeedLength, map3SeedMaxHits, map3SeedWindow, map3HPEnumeration )
122
123 def __main__():
124 #Parse Command Line
125 parser = optparse.OptionParser()
126 # Global options
127 parser.add_option( '--threads', dest='threads', help='The number of threads to use' )
128 parser.add_option( '--ref', dest='ref', help='The reference genome to use or index' )
129 parser.add_option( '--fastq', dest='fastq', help='The fastq file to use for the mapping' )
130 parser.add_option( '--output', dest='output', help='The file to save the output (SAM format)' )
131 parser.add_option( '--params', dest='params', help='Parameter setting to use (pre_set or full)' )
132 parser.add_option( '--fileSource', dest='fileSource', help='Whether to use a previously indexed reference sequence or one from history (indexed or history)' )
133 parser.add_option( '--mapall', dest='mapall', help='Flag indicating if mapall options are present')
134 parser.add_option( '--matchScore', dest='matchScore', help='The match score' )
135 parser.add_option( '--mismatchPenalty', dest='mismatchPenalty', help='Mismatch penalty' )
136 parser.add_option( '--gapOpenPenalty', dest='gapOpenPenalty', help='Gap open penalty' )
137 parser.add_option( '--gapExtensPenalty', dest='gapExtensPenalty', help='Gap extension penalty' )
138 parser.add_option( '--flowPenalty', dest='flowPenalty', help='Flow score penalty' )
139 parser.add_option( '--flowOrder', dest='flowOrder', help='Flow order' )
140 parser.add_option( '--bandWidth', dest='bandWidth', help='The band width' )
141 parser.add_option( '--globalMap', dest='globalMap', help='Map the full read (no soft-clipping)' )
142 parser.add_option( '--duplicateWindow', dest='duplicateWindow', help='Remove duplicate alignments from different algorithms within this bp window (-1 to disable)' )
143 parser.add_option( '--scoringThreshold', dest='scoringThreshold', help='The score threshold divided by the match score ' )
144 parser.add_option( '--queueSize', dest='queueSize', help='The queue size for the reads' )
145 parser.add_option( '--outputFilter', dest='outputFilter', help='The output filter (0 - unique best hits, 1 - random best hit, 2 - all best htis, 3 - all alignments)' )
146 parser.add_option( '--rgTag', dest='rgTag', help='The flag to specify RG tag(s) in the SAM header' )
147 parser.add_option( '--rgTagID', dest='rgTagID', default='', help='The RG ID to add to the SAM header' )
148 parser.add_option( '--rgTagCN', dest='rgTagCN', default='', help='The RG CN to add to the SAM header' )
149 parser.add_option( '--rgTagDS', dest='rgTagDS', default='', help='The RG DS to add to the SAM header' )
150 parser.add_option( '--rgTagDT', dest='rgTagDT', default='', help='The RG DT to add to the SAM header' )
151 parser.add_option( '--rgTagLB', dest='rgTagLB', default='', help='The RG LB to add to the SAM header' )
152 parser.add_option( '--rgTagPI', dest='rgTagPI', default='', help='The RG PI to add to the SAM header' )
153 parser.add_option( '--rgTagPL', dest='rgTagPL', default='', help='The RG PL to add to the SAM header' )
154 parser.add_option( '--rgTagPU', dest='rgTagPU', default='', help='The RG PU to add to the SAM header' )
155 parser.add_option( '--rgTagSM', dest='rgTagSM', default='', help='The RG SM to add to the SAM header' )
156 parser.add_option( '--filterIndependently', dest='filterIndependently', help='Apply the output filter for each algorithm independently' )
157 parser.add_option( '--suppressHeader', dest='suppressHeader', help='Suppress header' )
158 parser.add_option( '--dbkey', dest='dbkey', help='Dbkey for reference genome' )
159 parser.add_option( '--do_not_build_index', dest='do_not_build_index', action='store_true', help="Don't build index" )
160 # map 1 - stage 1
161 parser.add_option( '--map1', dest='map1', help='True if map1 should be run in the first stage' )
162 parser.add_option( '--map1SeedLength', dest='map1SeedLength', help='The k-mer length to seed CALs (-1 to disable)' )
163 parser.add_option( '--map1SeedMismatches', dest='map1SeedMismatches', help='The maximum number of mismatches in the seed ' )
164 parser.add_option( '--map1SecondarySeedLength', dest='map1SecondarySeedLength', help='The secondary seed length (-1 to disable)' )
165 parser.add_option( '--map1NumEdits', dest='map1NumEdits', help='The maximum number of edits or false-negative probability assuming the maximum error rate' )
166 parser.add_option( '--map1BaseError', dest='map1BaseError', help='The assumed per-base maximum error rate' )
167 parser.add_option( '--map1Mismatches', dest='map1Mismatches', help='The maximum number of or (read length) fraction of mismatches ' )
168 parser.add_option( '--map1GapOpens', dest='map1GapOpens', help='The maximum number of or (read length) fraction of indel starts ' )
169 parser.add_option( '--map1GapExtensions', dest='map1GapExtensions', help='The maximum number of or (read length) fraction of indel extensions ' )
170 parser.add_option( '--map1MaxCALsDeletion', dest='map1MaxCALsDeletion', help='The maximum number of CALs to extend a deletion ' )
171 parser.add_option( '--map1EndIndels', dest='map1EndIndels', help='Indels are not allowed within this number of bps from the end of the read ' )
172 parser.add_option( '--map1MaxOptimalCALs', dest='map1MaxOptimalCALs', help='Stop searching when INT optimal CALs have been found ' )
173 parser.add_option( '--map1MaxNodes', dest='map1MaxNodes', help='The maximum number of alignment nodes ' )
174 # map 2 - stage 1
175 parser.add_option( '--map2', dest='map2', help='True if map2 should be run in the first stage' )
176 parser.add_option( '--map2Coefficient', dest='map2Coefficient', help='The coefficient of length-threshold adjustment ' )
177 parser.add_option( '--map2SeedIntervalSize', dest='map2SeedIntervalSize', help='The maximum seeding interval size ' )
178 parser.add_option( '--map2ZBest', dest='map2ZBest', help='Keep the z-best nodes during prefix trie traversal' )
179 parser.add_option( '--map2ReverseTrigger', dest='map2ReverseTrigger', help='The # seeds to trigger reverse alignment ' )
180 # map 3 - stage 1
181 parser.add_option( '--map3', dest='map3', help='True if map3 should be run in the first stage' )
182 parser.add_option( '--map3SeedLength', dest='map3SeedLength', help='The k-mer length to seed CALs (-1 tunes to the genome size) ' )
183 parser.add_option( '--map3SeedMaxHits', dest='map3SeedMaxHits', help='The maximum number of hits returned by a seed ' )
184 parser.add_option( '--map3SeedWindow', dest='map3SeedWindow', help='The window of bases in which to group seeds ' )
185 parser.add_option( '--map3HPEnumeration', dest='map3HPEnumeration', help='The single homopolymer error difference for enumeration ' )
186 # map 1 - stage 2
187 parser.add_option( '--MAP1', dest='MAP1', help='True if map1 should be run in the second stage' )
188 parser.add_option( '--MAP1SeedLength', dest='MAP1SeedLength', help='The k-mer length to seed CALs (-1 to disable)' )
189 parser.add_option( '--MAP1SeedMismatches', dest='MAP1SeedMismatches', help='The maximum number of mismatches in the seed ' )
190 parser.add_option( '--MAP1SecondarySeedLength', dest='MAP1SecondarySeedLength', help='The secondary seed length (-1 to disable)' )
191 parser.add_option( '--MAP1NumEdits', dest='MAP1NumEdits', help='The maximum number of edits or false-negative probability assuming the maximum error rate' )
192 parser.add_option( '--MAP1BaseError', dest='MAP1BaseError', help='The assumed per-base maximum error rate' )
193 parser.add_option( '--MAP1Mismatches', dest='MAP1Mismatches', help='The maximum number of or (read length) fraction of mismatches ' )
194 parser.add_option( '--MAP1GapOpens', dest='MAP1GapOpens', help='The maximum number of or (read length) fraction of indel starts ' )
195 parser.add_option( '--MAP1GapExtensions', dest='MAP1GapExtensions', help='The maximum number of or (read length) fraction of indel extensions ' )
196 parser.add_option( '--MAP1MaxCALsDeletion', dest='MAP1MaxCALsDeletion', help='The maximum number of CALs to extend a deletion ' )
197 parser.add_option( '--MAP1EndIndels', dest='MAP1EndIndels', help='Indels are not allowed within this number of bps from the end of the read ' )
198 parser.add_option( '--MAP1MaxOptimalCALs', dest='MAP1MaxOptimalCALs', help='Stop searching when INT optimal CALs have been found ' )
199 parser.add_option( '--MAP1MaxNodes', dest='MAP1MaxNodes', help='The maximum number of alignment nodes ' )
200 # map 2 - stage 2
201 parser.add_option( '--MAP2', dest='MAP2', help='True if map2 should be run in the second stage' )
202 parser.add_option( '--MAP2Coefficient', dest='MAP2Coefficient', help='The coefficient of length-threshold adjustment ' )
203 parser.add_option( '--MAP2SeedIntervalSize', dest='MAP2SeedIntervalSize', help='The maximum seeding interval size ' )
204 parser.add_option( '--MAP2ZBest', dest='MAP2ZBest', help='Keep the z-best nodes during prefix trie traversal' )
205 parser.add_option( '--MAP2ReverseTrigger', dest='MAP2ReverseTrigger', help='The # seeds to trigger reverse alignment ' )
206 # map 3 - stage 2
207 parser.add_option( '--MAP3', dest='MAP3', help='True if map3 should be run in the second stage' )
208 parser.add_option( '--MAP3SeedLength', dest='MAP3SeedLength', help='The k-mer length to seed CALs (-1 tunes to the genome size) ' )
209 parser.add_option( '--MAP3SeedMaxHits', dest='MAP3SeedMaxHits', help='The maximum number of hits returned by a seed ' )
210 parser.add_option( '--MAP3SeedWindow', dest='MAP3SeedWindow', help='The window of bases in which to group seeds ' )
211 parser.add_option( '--MAP3HPEnumeration', dest='MAP3HPEnumeration', help='The single homopolymer error difference for enumeration ' )
212 # parse the options
213 (options, args) = parser.parse_args()
214
215 # output version # of tool
216 try:
217 tmp = tempfile.NamedTemporaryFile().name
218 tmp_stdout = open( tmp, 'wb' )
219 proc = subprocess.Popen( args='tmap --version 2>&1', shell=True, stdout=tmp_stdout )
220 tmp_stdout.close()
221 returncode = proc.wait()
222 stdout = None
223 for line in open( tmp_stdout.name, 'rb' ):
224 if line.lower().find( 'version' ) >= 0:
225 stdout = line.strip()
226 break
227 if stdout:
228 sys.stdout.write( 'TMAP %s\n' % stdout )
229 else:
230 raise Exception
231 except:
232 sys.stdout.write( 'Could not determine TMAP version\n' )
233
234 fastq = options.fastq
235
236 # make temp directory for placement of indices
237 tmp_index_dir = tempfile.mkdtemp()
238 tmp_dir = tempfile.mkdtemp()
239
240 # index if necessary
241 if options.fileSource == 'history' and not options.do_not_build_index:
242 ref_file = tempfile.NamedTemporaryFile( dir=tmp_index_dir )
243 ref_file_name = ref_file.name
244 ref_file.close()
245 os.symlink( options.ref, ref_file_name )
246 cmd1 = 'tmap index -f %s -v ' % ( ref_file_name )
247 try:
248 tmp = tempfile.NamedTemporaryFile( dir=tmp_index_dir ).name
249 tmp_stderr = open( tmp, 'wb' )
250 proc = subprocess.Popen( args=cmd1, shell=True, cwd=tmp_index_dir, stderr=tmp_stderr.fileno() )
251 returncode = proc.wait()
252 tmp_stderr.close()
253 # get stderr, allowing for case where it's very large
254 tmp_stderr = open( tmp, 'rb' )
255 stderr = ''
256 buffsize = 1048576
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 # clean up temp dirs
269 if os.path.exists( tmp_index_dir ):
270 shutil.rmtree( tmp_index_dir )
271 if os.path.exists( tmp_dir ):
272 shutil.rmtree( tmp_dir )
273 stop_err( 'Error indexing reference sequence. ' + str( e ) )
274 else:
275 ref_file_name = options.ref
276
277 # set up mapping and generate mapping command options
278 if options.params == 'pre_set':
279 mapall_options = '-n %s' % ( options.threads )
280 map1_options = 'map1'
281 map2_options = 'map2'
282 map3_options = 'map3'
283 MAP1_options = ''
284 MAP2_options = ''
285 MAP3_options = ''
286 else:
287 # mapall options
288 if options.rgTag == 'true':
289 rgTag = ''
290 if options.rgTagID != '':
291 rgTag += '-R "ID:' + options.rgTagID + '" '
292 if options.rgTagCN != '':
293 rgTag += '-R "CN:' + options.rgTagCN + '" '
294 if options.rgTagDS != '':
295 rgTag += '-R "DS:' + options.rgTagDS + '" '
296 if options.rgTagDT != '':
297 rgTag += '-R "DT:' + options.rgTagDT + '" '
298 if options.rgTagLB != '':
299 rgTag += '-R "LB:' + options.rgTagLB + '" '
300 if options.rgTagPI != '':
301 rgTag += '-R "PI:' + options.rgTagPI + '" '
302 if options.rgTagPL != '':
303 rgTag += '-R "PL:' + options.rgTagPL + '" '
304 if options.rgTagPU != '':
305 rgTag += '-R "PU:' + options.rgTagPU + '" '
306 if options.rgTagSM != '':
307 rgTag += '-R "SM:' + options.rgTagSM + '" '
308 rgTag.rstrip(' ')
309 else:
310 rgTag = ''
311 if None != options.flowOrder and '' != options.flowOrder:
312 flowOrder = '-x ' + options.flowOrder
313 else:
314 flowOrder = ''
315 if options.filterIndependently == 'true':
316 filterIndependently = '-I'
317 else:
318 filterIndependently = ''
319 if options.mapall == 'true':
320 mapall_options = '-A %s -M %s -O %s -E %s -X %s %s %s -W %s -T %s -q %s -n %s -a %s %s %s' % \
321 ( options.matchScore, options.mismatchPenalty, options.gapOpenPenalty, options.gapExtensPenalty,
322 options.flowPenalty, flowOrder,
323 options.globalMap, options.duplicateWindow, options.scoringThreshold, options.queueSize,
324 options.threads, options.outputFilter, rgTag, filterIndependently )
325 else:
326 mapall_options = ''
327
328 # map1 - stage one
329 if options.map1 == 'true':
330 map1_options = 'map1 %s' % \
331 map1_parse( options.map1SeedLength, options.map1SeedMismatches, \
332 options.map1SecondarySeedLength, options.map1NumEdits, \
333 options.map1BaseError, options.map1Mismatches, \
334 options.map1GapOpens, options.map1GapExtensions, options.map1MaxCALsDeletion, \
335 options.map1EndIndels, options.map1MaxOptimalCALs, options.map1MaxNodes )
336 else:
337 map1_options = ''
338 # map2 - stage one
339 if options.map2 == 'true':
340 map2_options = 'map2 %s' % \
341 map2_parse( options.map2Coefficient, options.map2SeedIntervalSize, \
342 options.map2ZBest, options.map2ReverseTrigger )
343 else:
344 map2_options = ''
345 # map3 - stage one
346 if options.map3 == 'true':
347 map3_options = 'map3 %s' % \
348 map3_parse( options.map3SeedLength, options.map3SeedMaxHits, \
349 options.map3SeedWindow, options.map3HPEnumeration )
350 else:
351 map3_options = ''
352 # map1 - stage two
353 if options.MAP1== 'true':
354 MAP1_options = 'MAP1 %s' % \
355 map1_parse( options.MAP1SeedLength, options.MAP1SeedMismatches, \
356 options.MAP1SecondarySeedLength, options.MAP1NumEdits, \
357 options.MAP1BaseError, options.MAP1Mismatches, \
358 options.MAP1GapOpens, options.MAP1GapExtensions, options.MAP1MaxCALsDeletion, \
359 options.MAP1EndIndels, options.MAP1MaxOptimalCALs, options.MAP1MaxNodes )
360 else:
361 MAP1_options = ''
362 # map2 - stage two
363 if options.MAP2 == 'true':
364 MAP2_options = 'MAP2 %s' % \
365 map2_parse( options.MAP2Coefficient, options.MAP2SeedIntervalSize, \
366 options.MAP2ZBest, options.MAP2ReverseTrigger )
367 else:
368 MAP2_options = ''
369 # map3 - stage two
370 if options.MAP3 == 'true':
371 MAP3_options = 'MAP3 %s' % \
372 MAP3_parse( options.MAP3SeedLength, options.MAP3SeedMaxHits, \
373 options.MAP3SeedWindow, options.MAP3HPEnumeration )
374 else:
375 MAP3_options = ''
376
377 #mapping_cmds
378 # prepare actual mapping and generate mapping commands
379 cmd2 = 'tmap mapall -f %s -r %s -F fastq %s -v %s %s %s %s %s %s > %s' % \
380 ( ref_file_name, fastq, mapall_options,
381 map1_options, map2_options, map3_options,
382 MAP1_options, MAP2_options, MAP3_options,
383 options.output )
384 # perform alignments
385 buffsize = 1048576
386 try:
387 # need to nest try-except in try-finally to handle 2.4
388 try:
389 # align
390 try:
391 tmp = tempfile.NamedTemporaryFile( dir=tmp_dir ).name
392 tmp_stderr = open( tmp, 'wb' )
393 proc = subprocess.Popen( args=cmd2, shell=True, cwd=tmp_dir, stderr=tmp_stderr.fileno() )
394 returncode = proc.wait()
395 tmp_stderr.close()
396 # get stderr, allowing for case where it's very large
397 tmp_stderr = open( tmp, 'rb' )
398 stderr = ''
399 try:
400 while True:
401 stderr += tmp_stderr.read( buffsize )
402 if not stderr or len( stderr ) % buffsize != 0:
403 break
404 except OverflowError:
405 pass
406 tmp_stderr.close()
407 if returncode != 0:
408 raise Exception, stderr
409 except Exception, e:
410 raise Exception, 'Error mapping sequence. ' + str( e )
411 # remove header if necessary
412 if options.suppressHeader == 'true':
413 tmp_out = tempfile.NamedTemporaryFile( dir=tmp_dir)
414 tmp_out_name = tmp_out.name
415 tmp_out.close()
416 try:
417 shutil.move( options.output, tmp_out_name )
418 except Exception, e:
419 raise Exception, 'Error moving output file before removing headers. ' + str( e )
420 fout = file( options.output, 'w' )
421 for line in file( tmp_out.name, 'r' ):
422 if not ( line.startswith( '@HD' ) or line.startswith( '@SQ' ) or line.startswith( '@RG' ) or line.startswith( '@PG' ) or line.startswith( '@CO' ) ):
423 fout.write( line )
424 fout.close()
425 # check that there are results in the output file
426 if os.path.getsize( options.output ) > 0:
427 sys.stdout.write( 'TMAP completed' )
428 else:
429 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.'
430 except Exception, e:
431 stop_err( 'The alignment failed.\n' + str( e ) )
432 finally:
433 # clean up temp dir
434 if os.path.exists( tmp_index_dir ):
435 shutil.rmtree( tmp_index_dir )
436 if os.path.exists( tmp_dir ):
437 shutil.rmtree( tmp_dir )
438
439 if __name__=="__main__": __main__()