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