Mercurial > repos > greg > gregs_test_repo
comparison galaxy_tmap_tool/tmap_wrapper.py @ 10:09951b7d29a4
Uploaded galaxy_tmap_tool.tar
| author | greg |
|---|---|
| date | Thu, 21 Jul 2011 09:54:55 -0400 |
| parents | |
| children |
comparison
equal
deleted
inserted
replaced
| 9:57045525eba1 | 10:09951b7d29a4 |
|---|---|
| 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__() |
