1 #include <stdlib.h>
2 #include <iostream>
3 #include <fstream>
4 #include <string>
5 #include <algorithm>
6 #include <cassert>
7 #include <stdexcept>
8 #include <getopt.h>
9 #include <vector>
10 #include <time.h>
11
12 #ifndef _WIN32
13 #include <dirent.h>
14 #include <signal.h>
15 #endif
16
17 #include "aligner.h"
18 #include "aligner_0mm.h"
19 #include "aligner_1mm.h"
20 #include "aligner_23mm.h"
21 #include "aligner_metrics.h"
22 #include "aligner_seed_mm.h"
23 #include "alphabet.h"
24 #include "assert_helpers.h"
25 #include "bitset.h"
26 #include "ds.h"
27 #include "ebwt.h"
28 #include "ebwt_search.h"
29 #include "endian_swap.h"
30 #include "formats.h"
31 #include "hit.h"
32 #include "pat.h"
33 #include "range_cache.h"
34 #include "sam.h"
35 #include "sequence_io.h"
36 #include "threading.h"
37 #include "tokenize.h"
38 #ifdef CHUD_PROFILING
39 #include <CHUD/CHUD.h>
40 #endif
41
42 static int FNAME_SIZE;
43 #if (__cplusplus >= 201103L)
44 #include <thread>
45 static std::atomic<int> thread_counter;
46 #else
47 static int thread_counter;
48 static MUTEX_T thread_counter_mutex;
49 #endif
50
51 using namespace std;
52
53 static EList<string> mates1; // mated reads (first mate)
54 static EList<string> mates2; // mated reads (second mate)
55 static EList<string> mates12; // mated reads (1st/2nd interleaved in 1 file)
56 static string adjustedEbwtFileBase;
57 static bool verbose; // be talkative
58 static bool startVerbose; // be talkative at startup
59 bool quiet; // print nothing but the alignments
60 static int sanityCheck; // enable expensive sanity checks
61 static int format; // default read format is FASTQ
62 static string origString; // reference text, or filename(s)
63 static int seed; // srandom() seed
64 static int timing; // whether to report basic timing data
65 static bool allHits; // for multihits, report just one
66 static bool rangeMode; // report BWT ranges instead of ref locs
67 static int showVersion; // just print version and quit?
68 static int ipause; // pause before maching?
69 static uint32_t qUpto; // max # of queries to read
70 static int trim5; // amount to trim from 5' end
71 static int trim3; // amount to trim from 3' end
72 static int reportOpps; // whether to report # of other mappings
73 static int offRate; // keep default offRate
74 static int isaRate; // keep default isaRate
75 static int mismatches; // allow 0 mismatches by default
76 static bool solexaQuals; // quality strings are solexa quals, not phred, and subtract 64 (not 33)
77 static bool phred64Quals; // quality chars are phred, but must subtract 64 (not 33)
78 static bool integerQuals; // quality strings are space-separated strings of integers, not ASCII
79 static int maqLike; // do maq-like searching
80 static int seedLen; // seed length (changed in Maq 0.6.4 from 24)
81 static int seedMms; // # mismatches allowed in seed (maq's -n)
82 static int qualThresh; // max qual-weighted hamming dist (maq's -e)
83 static int maxBtsBetter; // max # backtracks allowed in half-and-half mode
84 static int maxBts; // max # backtracks allowed in half-and-half mode
85 static int nthreads; // number of pthreads operating concurrently
86 static bool reorder; // reorder SAM output when running multi-threaded
87 static int thread_ceiling; // maximum number of threads user wants bowtie to use
88 static string thread_stealing_dir; // keep track of pids in this directory
89 static bool thread_stealing; // true iff thread stealing is in use
90 static output_types outType; // style of output
91 static bool noRefNames; // true -> print reference indexes; not names
92 static string dumpAlBase; // basename of same-format files to dump aligned reads to
93 static string dumpUnalBase; // basename of same-format files to dump unaligned reads to
94 static string dumpMaxBase; // basename of same-format files to dump reads with more than -m valid alignments to
95 static uint32_t khits; // number of hits per read; >1 is much slower
96 static uint32_t mhits; // don't report any hits if there are > mhits
97 static bool better; // true -> guarantee alignments from best possible stratum
98 static bool strata; // true -> don't stop at stratum boundaries
99 static int partitionSz; // output a partitioning key in first field
100 static int readsPerBatch; // # reads to read from input file at once
101 static size_t outBatchSz; // # alignments to write to output file at once
102 static bool noMaqRound; // true -> don't round quals to nearest 10 like maq
103 static bool fileParallel; // separate threads read separate input files in parallel
104 static bool useShmem; // use shared memory to hold the index
105 static bool useMm; // use memory-mapped files to hold the index
106 static bool mmSweep; // sweep through memory-mapped files immediately after mapping
107 static bool stateful; // use stateful aligners
108 static uint32_t prefetchWidth; // number of reads to process in parallel w/ --stateful
109 static uint32_t minInsert; // minimum insert size (Maq = 0, SOAP = 400)
110 static uint32_t maxInsert; // maximum insert size (Maq = 250, SOAP = 600)
111 static bool mate1fw; // -1 mate aligns in fw orientation on fw strand
112 static bool mate2fw; // -2 mate aligns in rc orientation on fw strand
113 static bool mateFwSet; // true -> user set --ff/--fr/--rf
114 static uint32_t mixedThresh; // threshold for when to switch to paired-end mixed mode (see aligner.h)
115 static uint32_t mixedAttemptLim; // number of attempts to make in "mixed mode" before giving up on orientation
116 static bool dontReconcileMates; // suppress pairwise all-versus-all way of resolving mates
117 static uint32_t cacheLimit; // ranges w/ size > limit will be cached
118 static uint32_t cacheSize; // # words per range cache
119 static int offBase; // offsets are 0-based by default, but configurable
120 static bool tryHard; // set very high maxBts, mixedAttemptLim
121 static uint32_t skipReads; // # reads/read pairs to skip
122 static bool nofw; // don't align fw orientation of read
123 static bool norc; // don't align rc orientation of read
124 static bool strandFix; // attempt to fix strand bias
125 static bool stats; // print performance stats
126 static int chunkPoolMegabytes; // max MB to dedicate to best-first search frames per thread
127 static int chunkSz; // size of single chunk disbursed by ChunkPool
128 static bool chunkVerbose; // have chunk allocator output status messages?
129 static bool useV1;
130 static bool reportSe;
131 static size_t fastaContLen;
132 static size_t fastaContFreq;
133 static bool hadoopOut; // print Hadoop status and summary messages
134 static bool fullRef;
135 static bool samNoQnameTrunc; // don't truncate QNAME field at first whitespace
136 static bool samNoHead; // don't print any header lines in SAM output
137 static bool samNoSQ; // don't print @SQ header lines
138 static string rgs; // SAM outputs for @RG header line
139 static Bitset suppressOuts(64); // output fields to suppress
140 static bool sampleMax; // whether to report a random alignment when maxed-out via -m/-M
141 static int defaultMapq; // default mapping quality to print in SAM mode
142 static bool printCost; // true -> print stratum and cost
143 bool showSeed;
144 static EList<string> qualities;
145 static EList<string> qualities1;
146 static EList<string> qualities2;
147 static string wrapper; // Type of wrapper script
148 bool gAllowMateContainment;
149 bool noUnal; // don't print unaligned reads
150 string ebwtFile; // read serialized Ebwt from this file
151 MUTEX_T gLock;
152
resetOptions()153 static void resetOptions() {
154 mates1.clear();
155 mates2.clear();
156 mates12.clear();
157 ebwtFile = "";
158 adjustedEbwtFileBase = "";
159 verbose = 0;
160 startVerbose = 0;
161 quiet = false;
162 sanityCheck = 0; // enable expensive sanity checks
163 format = FASTQ; // default read format is FASTQ
164 origString = ""; // reference text, or filename(s)
165 seed = 0; // srandom() seed
166 timing = 0; // whether to report basic timing data
167 allHits = false; // for multihits, report just one
168 rangeMode = false; // report BWT ranges instead of ref locs
169 showVersion = 0; // just print version and quit?
170 ipause = 0; // pause before maching?
171 qUpto = 0xffffffff; // max # of queries to read
172 trim5 = 0; // amount to trim from 5' end
173 trim3 = 0; // amount to trim from 3' end
174 reportOpps = 0; // whether to report # of other mappings
175 offRate = -1; // keep default offRate
176 isaRate = -1; // keep default isaRate
177 mismatches = 0; // allow 0 mismatches by default
178 solexaQuals = false; // quality strings are solexa quals, not phred, and subtract 64 (not 33)
179 phred64Quals = false; // quality chars are phred, but must subtract 64 (not 33)
180 integerQuals = false; // quality strings are space-separated strings of integers, not ASCII
181 maqLike = 1; // do maq-like searching
182 seedLen = 28; // seed length (changed in Maq 0.6.4 from 24)
183 seedMms = 2; // # mismatches allowed in seed (maq's -n)
184 qualThresh = 70; // max qual-weighted hamming dist (maq's -e)
185 maxBtsBetter = 125; // max # backtracks allowed in half-and-half mode
186 maxBts = 800; // max # backtracks allowed in half-and-half mode
187 nthreads = 1; // number of pthreads operating concurrently
188 reorder = false; // reorder SAM output
189 thread_ceiling = 0; // max # threads user asked for
190 thread_stealing_dir = ""; // keep track of pids in this directory
191 thread_stealing = false; // true iff thread stealing is in use
192 FNAME_SIZE = 200;
193 outType = OUTPUT_FULL; // style of output
194 noRefNames = false; // true -> print reference indexes; not names
195 dumpAlBase = ""; // basename of same-format files to dump aligned reads to
196 dumpUnalBase = ""; // basename of same-format files to dump unaligned reads to
197 dumpMaxBase = ""; // basename of same-format files to dump reads with more than -m valid alignments to
198 khits = 1; // number of hits per read; >1 is much slower
199 mhits = 0xffffffff; // don't report any hits if there are > mhits
200 better = false; // true -> guarantee alignments from best possible stratum
201 strata = false; // true -> don't stop at stratum boundaries
202 partitionSz = 0; // output a partitioning key in first field
203 readsPerBatch = 16; // # reads to read from input file at once
204 outBatchSz = 16; // # alignments to wrote to output file at once
205 noMaqRound = false; // true -> don't round quals to nearest 10 like maq
206 fileParallel = false; // separate threads read separate input files in parallel
207 useShmem = false; // use shared memory to hold the index
208 useMm = false; // use memory-mapped files to hold the index
209 mmSweep = false; // sweep through memory-mapped files immediately after mapping
210 stateful = false; // use stateful aligners
211 prefetchWidth = 1; // number of reads to process in parallel w/ --stateful
212 minInsert = 0; // minimum insert size (Maq = 0, SOAP = 400)
213 maxInsert = 250; // maximum insert size (Maq = 250, SOAP = 600)
214 mate1fw = true; // -1 mate aligns in fw orientation on fw strand
215 mate2fw = false; // -2 mate aligns in rc orientation on fw strand
216 mateFwSet = false; // true -> user set mate1fw/mate2fw with --ff/--fr/--rf
217 mixedThresh = 4; // threshold for when to switch to paired-end mixed mode (see aligner.h)
218 mixedAttemptLim = 100; // number of attempts to make in "mixed mode" before giving up on orientation
219 dontReconcileMates = true; // suppress pairwise all-versus-all way of resolving mates
220 cacheLimit = 5; // ranges w/ size > limit will be cached
221 cacheSize = 0; // # words per range cache
222 offBase = 0; // offsets are 0-based by default, but configurable
223 tryHard = false; // set very high maxBts, mixedAttemptLim
224 skipReads = 0; // # reads/read pairs to skip
225 nofw = false; // don't align fw orientation of read
226 norc = false; // don't align rc orientation of read
227 strandFix = true; // attempt to fix strand bias
228 stats = false; // print performance stats
229 chunkPoolMegabytes = 64; // max MB to dedicate to best-first search frames per thread
230 chunkSz = 256; // size of single chunk disbursed by ChunkPool (in KB)
231 chunkVerbose = false; // have chunk allocator output status messages?
232 useV1 = true;
233 reportSe = false;
234 fastaContLen = 0;
235 fastaContFreq = 0;
236 hadoopOut = false; // print Hadoop status and summary messages
237 fullRef = false; // print entire reference name instead of just up to 1st space
238 samNoQnameTrunc = false; // don't truncate at first whitespace?
239 samNoHead = false; // don't print any header lines in SAM output
240 samNoSQ = false; // don't print @SQ header lines
241 rgs = ""; // SAM outputs for @RG header line
242 suppressOuts.clear(); // output fields to suppress
243 sampleMax = false;
244 defaultMapq = 255;
245 printCost = false; // true -> print cost and stratum
246 showSeed = false; // true -> print per-read pseudo-random seed
247 qualities.clear();
248 qualities1.clear();
249 qualities2.clear();
250 wrapper.clear();
251 gAllowMateContainment = false; // true -> alignments where one mate lies inside the other are valid
252 noUnal = false; // true -> do not report unaligned reads
253 }
254
255 // mating constraints
256
257 static const char *short_options = "fF:qbzhcu:rv:s:at3:5:o:e:n:l:w:p:k:m:M:1:2:I:X:x:z:B:ySCQ:";
258
259 enum {
260 ARG_ORIG = 256,
261 ARG_SEED,
262 ARG_RANGE,
263 ARG_SOLEXA_QUALS,
264 ARG_MAXBTS,
265 ARG_VERBOSE,
266 ARG_STARTVERBOSE,
267 ARG_QUIET,
268 ARG_FAST,
269 ARG_AL,
270 ARG_UN,
271 ARG_MAXDUMP,
272 ARG_REFIDX,
273 ARG_SANITY,
274 ARG_OLDBEST,
275 ARG_BETTER,
276 ARG_BEST,
277 ARG_ISARATE,
278 ARG_PARTITION,
279 ARG_READS_PER_BATCH,
280 ARG_integerQuals,
281 ARG_NOMAQROUND,
282 ARG_FILEPAR,
283 ARG_SHMEM,
284 ARG_MM,
285 ARG_MMSWEEP,
286 ARG_STATEFUL,
287 ARG_PREFETCH_WIDTH,
288 ARG_FF,
289 ARG_FR,
290 ARG_RF,
291 ARG_MIXED_ATTEMPTS,
292 ARG_NO_RECONCILE,
293 ARG_CACHE_LIM,
294 ARG_CACHE_SZ,
295 ARG_NO_FW,
296 ARG_NO_RC,
297 ARG_SKIP,
298 ARG_STRAND_FIX,
299 ARG_STATS,
300 ARG_ONETWO,
301 ARG_PHRED64,
302 ARG_PHRED33,
303 ARG_CHUNKMBS,
304 ARG_CHUNKSZ,
305 ARG_CHUNKVERBOSE,
306 ARG_STRATA,
307 ARG_PEV2,
308 ARG_REPORTSE,
309 ARG_HADOOPOUT,
310 ARG_FUZZY,
311 ARG_FULLREF,
312 ARG_USAGE,
313 ARG_SAM_NO_QNAME_TRUNC,
314 ARG_SAM_NOHEAD,
315 ARG_SAM_NOSQ,
316 ARG_SAM_RG,
317 ARG_SUPPRESS_FIELDS,
318 ARG_DEFAULT_MAPQ,
319 ARG_COST,
320 ARG_SHOWSEED,
321 ARG_QUALS1,
322 ARG_QUALS2,
323 ARG_ALLOW_CONTAIN,
324 ARG_WRAPPER,
325 ARG_INTERLEAVED_FASTQ,
326 ARG_SAM_NO_UNAL,
327 ARG_THREAD_CEILING,
328 ARG_THREAD_PIDDIR,
329 ARG_REORDER_SAM,
330 };
331
332 static struct option long_options[] = {
333 {(char*)"verbose", no_argument, 0, ARG_VERBOSE},
334 {(char*)"startverbose", no_argument, 0, ARG_STARTVERBOSE},
335 {(char*)"quiet", no_argument, 0, ARG_QUIET},
336 {(char*)"sanity", no_argument, 0, ARG_SANITY},
337 {(char*)"pause", no_argument, &ipause, 1},
338 {(char*)"orig", required_argument, 0, ARG_ORIG},
339 {(char*)"all", no_argument, 0, 'a'},
340 {(char*)"solexa-quals", no_argument, 0, ARG_SOLEXA_QUALS},
341 {(char*)"integer-quals", no_argument, 0, ARG_integerQuals},
342 {(char*)"time", no_argument, 0, 't'},
343 {(char*)"trim3", required_argument, 0, '3'},
344 {(char*)"trim5", required_argument, 0, '5'},
345 {(char*)"seed", required_argument, 0, ARG_SEED},
346 {(char*)"qupto", required_argument, 0, 'u'},
347 {(char*)"al", required_argument, 0, ARG_AL},
348 {(char*)"un", required_argument, 0, ARG_UN},
349 {(char*)"max", required_argument, 0, ARG_MAXDUMP},
350 {(char*)"offrate", required_argument, 0, 'o'},
351 {(char*)"isarate", required_argument, 0, ARG_ISARATE},
352 {(char*)"reportopps", no_argument, &reportOpps, 1},
353 {(char*)"version", no_argument, &showVersion, 1},
354 {(char*)"reads-per-batch", required_argument, 0, ARG_READS_PER_BATCH},
355 {(char*)"maqerr", required_argument, 0, 'e'},
356 {(char*)"seedlen", required_argument, 0, 'l'},
357 {(char*)"seedmms", required_argument, 0, 'n'},
358 {(char*)"filepar", no_argument, 0, ARG_FILEPAR},
359 {(char*)"help", no_argument, 0, 'h'},
360 {(char*)"threads", required_argument, 0, 'p'},
361 {(char*)"khits", required_argument, 0, 'k'},
362 {(char*)"mhits", required_argument, 0, 'm'},
363 {(char*)"minins", required_argument, 0, 'I'},
364 {(char*)"maxins", required_argument, 0, 'X'},
365 {(char*)"quals", required_argument, 0, 'Q'},
366 {(char*)"Q1", required_argument, 0, ARG_QUALS1},
367 {(char*)"Q2", required_argument, 0, ARG_QUALS2},
368 {(char*)"best", no_argument, 0, ARG_BEST},
369 {(char*)"better", no_argument, 0, ARG_BETTER},
370 {(char*)"oldbest", no_argument, 0, ARG_OLDBEST},
371 {(char*)"strata", no_argument, 0, ARG_STRATA},
372 {(char*)"nomaqround", no_argument, 0, ARG_NOMAQROUND},
373 {(char*)"refidx", no_argument, 0, ARG_REFIDX},
374 {(char*)"range", no_argument, 0, ARG_RANGE},
375 {(char*)"maxbts", required_argument, 0, ARG_MAXBTS},
376 {(char*)"phased", no_argument, 0, 'z'},
377 {(char*)"partition", required_argument, 0, ARG_PARTITION},
378 {(char*)"stateful", no_argument, 0, ARG_STATEFUL},
379 {(char*)"prewidth", required_argument, 0, ARG_PREFETCH_WIDTH},
380 {(char*)"ff", no_argument, 0, ARG_FF},
381 {(char*)"fr", no_argument, 0, ARG_FR},
382 {(char*)"rf", no_argument, 0, ARG_RF},
383 {(char*)"mixthresh", required_argument, 0, 'x'},
384 {(char*)"pairtries", required_argument, 0, ARG_MIXED_ATTEMPTS},
385 {(char*)"noreconcile", no_argument, 0, ARG_NO_RECONCILE},
386 {(char*)"cachelim", required_argument, 0, ARG_CACHE_LIM},
387 {(char*)"cachesz", required_argument, 0, ARG_CACHE_SZ},
388 {(char*)"nofw", no_argument, 0, ARG_NO_FW},
389 {(char*)"norc", no_argument, 0, ARG_NO_RC},
390 {(char*)"offbase", required_argument, 0, 'B'},
391 {(char*)"tryhard", no_argument, 0, 'y'},
392 {(char*)"skip", required_argument, 0, 's'},
393 {(char*)"strandfix", no_argument, 0, ARG_STRAND_FIX},
394 {(char*)"stats", no_argument, 0, ARG_STATS},
395 {(char*)"12", required_argument, 0, ARG_ONETWO},
396 {(char*)"phred33-quals", no_argument, 0, ARG_PHRED33},
397 {(char*)"phred64-quals", no_argument, 0, ARG_PHRED64},
398 {(char*)"solexa1.3-quals", no_argument, 0, ARG_PHRED64},
399 {(char*)"chunkmbs", required_argument, 0, ARG_CHUNKMBS},
400 {(char*)"chunksz", required_argument, 0, ARG_CHUNKSZ},
401 {(char*)"chunkverbose", no_argument, 0, ARG_CHUNKVERBOSE},
402 {(char*)"mm", no_argument, 0, ARG_MM},
403 {(char*)"shmem", no_argument, 0, ARG_SHMEM},
404 {(char*)"mmsweep", no_argument, 0, ARG_MMSWEEP},
405 {(char*)"pev2", no_argument, 0, ARG_PEV2},
406 {(char*)"reportse", no_argument, 0, ARG_REPORTSE},
407 {(char*)"hadoopout", no_argument, 0, ARG_HADOOPOUT},
408 {(char*)"fullref", no_argument, 0, ARG_FULLREF},
409 {(char*)"usage", no_argument, 0, ARG_USAGE},
410 {(char*)"sam", no_argument, 0, 'S'},
411 {(char*)"sam-no-qname-trunc", no_argument, 0, ARG_SAM_NO_QNAME_TRUNC},
412 {(char*)"sam-nohead", no_argument, 0, ARG_SAM_NOHEAD},
413 {(char*)"sam-nosq", no_argument, 0, ARG_SAM_NOSQ},
414 {(char*)"sam-noSQ", no_argument, 0, ARG_SAM_NOSQ},
415 {(char*)"sam-RG", required_argument, 0, ARG_SAM_RG},
416 {(char*)"suppress", required_argument, 0, ARG_SUPPRESS_FIELDS},
417 {(char*)"mapq", required_argument, 0, ARG_DEFAULT_MAPQ},
418 {(char*)"cost", no_argument, 0, ARG_COST},
419 {(char*)"showseed", no_argument, 0, ARG_SHOWSEED},
420 {(char*)"allow-contain",no_argument, 0, ARG_ALLOW_CONTAIN},
421 {(char*)"wrapper", required_argument, 0, ARG_WRAPPER},
422 {(char*)"interleaved", required_argument, 0, ARG_INTERLEAVED_FASTQ},
423 {(char*)"no-unal", no_argument, 0, ARG_SAM_NO_UNAL},
424 {(char*)"thread-ceiling",required_argument, 0, ARG_THREAD_CEILING},
425 {(char*)"thread-piddir", required_argument, 0, ARG_THREAD_PIDDIR},
426 {(char*)"reorder", no_argument, 0, ARG_REORDER_SAM},
427 {(char*)0, 0, 0, 0} // terminator
428 };
429
430 /**
431 * Print a summary usage message to the provided output stream.
432 */
printUsage(ostream & out)433 static void printUsage(ostream& out) {
434 #ifdef BOWTIE_64BIT_INDEX
435 string tool_name = "bowtie-align-l";
436 #else
437 string tool_name = "bowtie-align-s";
438 #endif
439 if(wrapper == "basic-0") {
440 tool_name = "bowtie";
441 }
442
443 out << "Usage: " << endl
444 << tool_name << " [options]* -x <ebwt> {-1 <m1> -2 <m2> | --12 <r> | --interleaved <i> | <s>} [<hit>]" << endl
445 << endl
446 << " <ebwt> Index filename prefix (minus trailing .X." + gEbwt_ext + ")." << endl
447 << " <m1> Comma-separated list of files containing upstream mates (or the" << endl
448 << " sequences themselves, if -c is set) paired with mates in <m2>" << endl
449 << " <m2> Comma-separated list of files containing downstream mates (or the" << endl
450 << " sequences themselves if -c is set) paired with mates in <m1>" << endl
451 << " <r> Comma-separated list of files containing Crossbow-style reads. Can be" << endl
452 << " a mixture of paired and unpaired. Specify \"-\" for stdin." << endl
453 << " <i> Files with interleaved paired-end FASTQ reads." << endl
454 << " <s> Comma-separated list of files containing unpaired reads, or the" << endl
455 << " sequences themselves, if -c is set. Specify \"-\" for stdin." << endl
456 << " <hit> File to write hits to (default: stdout)" << endl
457 << "Input:" << endl
458 << " -q query input files are FASTQ .fq/.fastq (default)" << endl
459 << " -f query input files are (multi-)FASTA .fa/.mfa" << endl
460 << " -F k:<int>,i:<int> query input files are continuous FASTA where reads" << endl
461 << " are substrings (k-mers) extracted from a FASTA file <s>" << endl
462 << " and aligned at offsets 1, 1+i, 1+2i ... end of reference" << endl
463 << " -r query input files are raw one-sequence-per-line" << endl
464 << " -c query sequences given on cmd line (as <mates>, <singles>)" << endl
465 << " -Q/--quals <file> QV file(s) corresponding to CSFASTA inputs; use with -f -C" << endl
466 << " --Q1/--Q2 <file> same as -Q, but for mate files 1 and 2 respectively" << endl
467 << " -s/--skip <int> skip the first <int> reads/pairs in the input" << endl
468 << " -u/--qupto <int> stop after first <int> reads/pairs (excl. skipped reads)" << endl
469 << " -5/--trim5 <int> trim <int> bases from 5' (left) end of reads" << endl
470 << " -3/--trim3 <int> trim <int> bases from 3' (right) end of reads" << endl
471 << " --phred33-quals input quals are Phred+33 (default)" << endl
472 << " --phred64-quals input quals are Phred+64 (same as --solexa1.3-quals)" << endl
473 << " --solexa-quals input quals are from GA Pipeline ver. < 1.3" << endl
474 << " --solexa1.3-quals input quals are from GA Pipeline ver. >= 1.3" << endl
475 << " --integer-quals qualities are given as space-separated integers (not ASCII)" << endl;
476 if(wrapper == "basic-0") {
477 out << " --large-index force usage of a 'large' index, even if a small one is present" << endl;
478 }
479 out << "Alignment:" << endl
480 << " -v <int> report end-to-end hits w/ <=v mismatches; ignore qualities" << endl
481 << " or" << endl
482 << " -n/--seedmms <int> max mismatches in seed (can be 0-3, default: -n 2)" << endl
483 << " -e/--maqerr <int> max sum of mismatch quals across alignment for -n (def: 70)" << endl
484 << " -l/--seedlen <int> seed length for -n (default: 28)" << endl
485 << " --nomaqround disable Maq-like quality rounding for -n (nearest 10 <= 30)" << endl
486 << " -I/--minins <int> minimum insert size for paired-end alignment (default: 0)" << endl
487 << " -X/--maxins <int> maximum insert size for paired-end alignment (default: 250)" << endl
488 << " --fr/--rf/--ff -1, -2 mates align fw/rev, rev/fw, fw/fw (default: --fr)" << endl
489 << " --nofw/--norc do not align to forward/reverse-complement reference strand" << endl
490 << " --maxbts <int> max # backtracks for -n 2/3 (default: 125, 800 for --best)" << endl
491 << " --pairtries <int> max # attempts to find mate for anchor hit (default: 100)" << endl
492 << " -y/--tryhard try hard to find valid alignments, at the expense of speed" << endl
493 << " --chunkmbs <int> max megabytes of RAM for best-first search frames (def: 64)" << endl
494 << " --reads-per-batch # of reads to read from input file at once (default: 16)" << endl
495 << "Reporting:" << endl
496 << " -k <int> report up to <int> good alignments per read (default: 1)" << endl
497 << " -a/--all report all alignments per read (much slower than low -k)" << endl
498 << " -m <int> suppress all alignments if > <int> exist (def: no limit)" << endl
499 << " -M <int> like -m, but reports 1 random hit (MAPQ=0); requires --best" << endl
500 << " --best hits guaranteed best stratum; ties broken by quality" << endl
501 << " --strata hits in sub-optimal strata aren't reported (requires --best)" << endl
502 << "Output:" << endl
503 << " -t/--time print wall-clock time taken by search phases" << endl
504 << " -B/--offbase <int> leftmost ref offset = <int> in bowtie output (default: 0)" << endl
505 << " --quiet print nothing but the alignments" << endl
506 << " --refidx refer to ref. seqs by 0-based index rather than name" << endl
507 << " --al <fname> write aligned reads/pairs to file(s) <fname>" << endl
508 << " --un <fname> write unaligned reads/pairs to file(s) <fname>" << endl
509 << " --no-unal suppress SAM records for unaligned reads" << endl
510 << " --max <fname> write reads/pairs over -m limit to file(s) <fname>" << endl
511 << " --suppress <cols> suppresses given columns (comma-delim'ed) in default output" << endl
512 << " --fullref write entire ref name (default: only up to 1st space)" << endl
513 << "SAM:" << endl
514 << " -S/--sam write hits in SAM format" << endl
515 << " --mapq <int> default mapping quality (MAPQ) to print for SAM alignments" << endl
516 << " --sam-nohead supppress header lines (starting with @) for SAM output" << endl
517 << " --sam-nosq supppress @SQ header lines for SAM output" << endl
518 << " --sam-RG <text> add <text> (usually \"lab=value\") to @RG line of SAM header" << endl
519 << "Performance:" << endl
520 << " -o/--offrate <int> override offrate of index; must be >= index's offrate" << endl
521 << " -p/--threads <int> number of alignment threads to launch (default: 1)" << endl
522 #ifdef BOWTIE_MM
523 << " --mm use memory-mapped I/O for index; many 'bowtie's can share" << endl
524 #endif
525 #ifdef BOWTIE_SHARED_MEM
526 << " --shmem use shared mem for index; many 'bowtie's can share" << endl
527 #endif
528 << "Other:" << endl
529 << " --seed <int> seed for random number generator" << endl
530 << " --verbose verbose output (for debugging)" << endl
531 << " --version print version information and quit" << endl
532 << " -h/--help print this usage message" << endl
533 ;
534 if(wrapper.empty()) {
535 cerr << endl
536 << "*** Warning ***" << endl
537 << tool_name << " was run directly. It is recommended that you run the wrapper script 'bowtie' instead." << endl
538 << endl;
539 }
540 }
541
542 /**
543 * Parse an int out of optarg and enforce that it be at least 'lower';
544 * if it is less than 'lower', than output the given error message and
545 * exit with an error and a usage message.
546 */
parseInt(int lower,int upper,const char * errmsg,const char * arg)547 static int parseInt(int lower, int upper, const char *errmsg, const char *arg) {
548 long l;
549 char *endPtr= NULL;
550 l = strtol(arg, &endPtr, 10);
551 if (endPtr != NULL) {
552 if (l < lower || l > upper) {
553 cerr << errmsg << endl;
554 printUsage(cerr);
555 throw 1;
556 }
557 return (int32_t)l;
558 }
559 cerr << errmsg << endl;
560 printUsage(cerr);
561 throw 1;
562 return -1;
563 }
564
565 /**
566 * Parse from optarg by default.
567 */
parseInt(int lower,const char * errmsg)568 static int parseInt(int lower, const char *errmsg) {
569 return parseInt(lower, INT_MAX, errmsg, optarg);
570 }
571
572 /**
573 * Upper is INT_MAX by default.
574 */
parseInt(int lower,const char * errmsg,const char * arg)575 static int parseInt(int lower, const char *errmsg, const char *arg) {
576 return parseInt(lower, INT_MAX, errmsg, arg);
577 }
578
579 /**
580 * Upper is INT_MAX, parse from optarg by default.
581 */
parseInt(int lower,int upper,const char * errmsg)582 static int parseInt(int lower, int upper, const char *errmsg) {
583 return parseInt(lower, upper, errmsg, optarg);
584 }
585
586 /**
587 * Parse a T string 'str'.
588 */
589 template<typename T>
parse(const char * s)590 T parse(const char *s) {
591 T tmp;
592 stringstream ss(s);
593 ss >> tmp;
594 return tmp;
595 }
596
597 /**
598 * Parse a pair of Ts from a string, 'str', delimited with 'delim'.
599 */
600 template<typename T>
parsePair(const char * str,char delim)601 pair<T, T> parsePair(const char *str, char delim) {
602 string s(str);
603 EList<string> ss;
604 tokenize(s, delim, ss);
605 pair<T, T> ret;
606 ret.first = parse<T>(ss[0].c_str());
607 ret.second = parse<T>(ss[1].c_str());
608 return ret;
609 }
610
611 /**
612 * Read command-line arguments
613 */
parseOptions(int argc,const char ** argv)614 static void parseOptions(int argc, const char **argv) {
615 int option_index = 0;
616 int next_option;
617 if(startVerbose) { cerr << "Parsing options: "; logTime(cerr, true); }
618 do {
619 next_option = getopt_long(
620 argc, const_cast<char**>(argv),
621 short_options, long_options, &option_index);
622 switch (next_option) {
623 case ARG_WRAPPER: wrapper = optarg; break;
624 case '1': tokenize(optarg, ",", mates1); break;
625 case '2': tokenize(optarg, ",", mates2); break;
626 case ARG_ONETWO: tokenize(optarg, ",", mates12); format = TAB_MATE; break;
627 case ARG_INTERLEAVED_FASTQ: tokenize(optarg, ",", mates12); format = INTERLEAVED; break;
628 case 'f': format = FASTA; break;
629 case 'F': {
630 format = FASTA_CONT;
631 pair<size_t, size_t> p = parsePair<size_t>(optarg, ',');
632 fastaContLen = p.first;
633 fastaContFreq = p.second;
634 break;
635 }
636 case 'q': format = FASTQ; break;
637 case 'r': format = RAW; break;
638 case 'c': format = CMDLINE; break;
639 case 'I':
640 minInsert = (uint32_t)parseInt(0, "-I arg must be positive");
641 break;
642 case 'X':
643 maxInsert = (uint32_t)parseInt(1, "-X arg must be at least 1");
644 break;
645 case 'x':
646 ebwtFile = optarg;
647 break;
648 case 's':
649 skipReads = (uint32_t)parseInt(0, "-s arg must be positive");
650 break;
651 case ARG_FF: mate1fw = true; mate2fw = true; mateFwSet = true; break;
652 case ARG_RF: mate1fw = false; mate2fw = true; mateFwSet = true; break;
653 case ARG_FR: mate1fw = true; mate2fw = false; mateFwSet = true; break;
654 case ARG_RANGE: rangeMode = true; break;
655 case 'S': outType = OUTPUT_SAM; break;
656 case ARG_SHMEM: useShmem = true; break;
657 case ARG_SHOWSEED: showSeed = true; break;
658 case ARG_ALLOW_CONTAIN: gAllowMateContainment = true; break;
659 case ARG_SUPPRESS_FIELDS: {
660 EList<string> supp;
661 tokenize(optarg, ",", supp);
662 for(size_t i = 0; i < supp.size(); i++) {
663 int ii = parseInt(1, "--suppress arg must be at least 1", supp[i].c_str());
664 suppressOuts.set(ii-1);
665 }
666 break;
667 }
668 case ARG_MM: {
669 #ifdef BOWTIE_MM
670 useMm = true;
671 break;
672 #else
673 cerr << "Memory-mapped I/O mode is disabled because bowtie was not compiled with" << endl
674 << "BOWTIE_MM defined. Memory-mapped I/O is not supported under Windows. If you" << endl
675 << "would like to use memory-mapped I/O on a platform that supports it, please" << endl
676 << "refrain from specifying BOWTIE_MM=0 when compiling Bowtie." << endl;
677 throw 1;
678 #endif
679 }
680 case ARG_MMSWEEP: mmSweep = true; break;
681 case ARG_HADOOPOUT: hadoopOut = true; break;
682 case ARG_AL: dumpAlBase = optarg; break;
683 case ARG_UN: dumpUnalBase = optarg; break;
684 case ARG_MAXDUMP: dumpMaxBase = optarg; break;
685 case ARG_SOLEXA_QUALS: solexaQuals = true; break;
686 case ARG_integerQuals: integerQuals = true; break;
687 case ARG_PHRED64: phred64Quals = true; break;
688 case ARG_PHRED33: solexaQuals = false; phred64Quals = false; break;
689 case ARG_NOMAQROUND: noMaqRound = true; break;
690 case ARG_REFIDX: noRefNames = true; break;
691 case ARG_STATEFUL: stateful = true; break;
692 case ARG_REPORTSE: reportSe = true; break;
693 case ARG_FULLREF: fullRef = true; break;
694 case ARG_PREFETCH_WIDTH:
695 prefetchWidth = parseInt(1, "--prewidth must be at least 1");
696 break;
697 case 'B':
698 offBase = parseInt(-999999, "-B/--offbase cannot be a large negative number");
699 break;
700 case ARG_SEED:
701 seed = parseInt(0, "--seed arg must be at least 0");
702 break;
703 case 'u':
704 qUpto = (uint32_t)parseInt(1, "-u/--qupto arg must be at least 1");
705 break;
706 case 'k':
707 khits = (uint32_t)parseInt(1, "-k arg must be at least 1");
708 break;
709 case 'Q':
710 tokenize(optarg, ",", qualities);
711 integerQuals = true;
712 break;
713 case ARG_QUALS1:
714 tokenize(optarg, ",", qualities1);
715 integerQuals = true;
716 break;
717 case ARG_QUALS2:
718 tokenize(optarg, ",", qualities2);
719 integerQuals = true;
720 break;
721 case 'M':
722 sampleMax = true;
723 case 'm':
724 mhits = (uint32_t)parseInt(1, "-m arg must be at least 1");
725 break;
726 case 'z':
727 mixedThresh = (uint32_t)parseInt(0, "-x arg must be at least 0");
728 break;
729 case ARG_MIXED_ATTEMPTS:
730 mixedAttemptLim = (uint32_t)parseInt(1, "--pairtries arg must be at least 1");
731 break;
732 case ARG_CACHE_LIM:
733 cacheLimit = (uint32_t)parseInt(1, "--cachelim arg must be at least 1");
734 break;
735 case ARG_CACHE_SZ:
736 cacheSize = (uint32_t)parseInt(1, "--cachesz arg must be at least 1");
737 cacheSize *= (1024 * 1024); // convert from MB to B
738 break;
739 case ARG_NO_RECONCILE:
740 dontReconcileMates = true;
741 break;
742 case 'p':
743 nthreads = parseInt(1, "-p/--threads arg must be at least 1");
744 break;
745 case ARG_THREAD_CEILING:
746 thread_ceiling = parseInt(0, "--thread-ceiling must be at least 0");
747 break;
748 case ARG_THREAD_PIDDIR:
749 thread_stealing_dir = optarg;
750 break;
751 case ARG_REORDER_SAM:
752 reorder = true;
753 break;
754 case ARG_FILEPAR:
755 fileParallel = true;
756 break;
757 case 'v':
758 maqLike = 0;
759 mismatches = parseInt(0, 3, "-v arg must be at least 0 and at most 3");
760 break;
761 case '3': trim3 = parseInt(0, "-3/--trim3 arg must be at least 0"); break;
762 case '5': trim5 = parseInt(0, "-5/--trim5 arg must be at least 0"); break;
763 case 'o': offRate = parseInt(1, "-o/--offrate arg must be at least 1"); break;
764 case ARG_ISARATE: isaRate = parseInt(0, "--isarate arg must be at least 0"); break;
765 case 'e': qualThresh = parseInt(1, "-e/--err arg must be at least 1"); break;
766 case 'n': seedMms = parseInt(0, 3, "-n/--seedmms arg must be at least 0 and at most 3"); maqLike = 1; break;
767 case 'l': seedLen = parseInt(5, "-l/--seedlen arg must be at least 5"); break;
768 case 'h': printUsage(cout); throw 0; break;
769 case ARG_USAGE: printUsage(cout); throw 0; break;
770 case 'a': allHits = true; break;
771 case 'y': tryHard = true; break;
772 case ARG_CHUNKMBS: chunkPoolMegabytes = parseInt(1, "--chunkmbs arg must be at least 1"); break;
773 case ARG_CHUNKSZ: chunkSz = parseInt(1, "--chunksz arg must be at least 1"); break;
774 case ARG_CHUNKVERBOSE: chunkVerbose = true; break;
775 case ARG_BETTER: stateful = true; better = true; break;
776 case ARG_BEST: stateful = true; useV1 = false; break;
777 case ARG_STRATA: strata = true; break;
778 case ARG_VERBOSE: verbose = true; break;
779 case ARG_STARTVERBOSE: startVerbose = true; break;
780 case ARG_QUIET: quiet = true; break;
781 case ARG_SANITY: sanityCheck = true; break;
782 case 't': timing = true; break;
783 case ARG_NO_FW: nofw = true; break;
784 case ARG_NO_RC: norc = true; break;
785 case ARG_STATS: stats = true; break;
786 case ARG_PEV2: useV1 = false; break;
787 case ARG_SAM_NO_QNAME_TRUNC: samNoQnameTrunc = true; break;
788 case ARG_SAM_NOHEAD: samNoHead = true; break;
789 case ARG_SAM_NOSQ: samNoSQ = true; break;
790 case ARG_SAM_NO_UNAL: noUnal = true; break;
791 case ARG_SAM_RG: {
792 if(!rgs.empty()) rgs += '\t';
793 rgs += optarg;
794 break;
795 }
796 case ARG_COST: printCost = true; break;
797 case ARG_DEFAULT_MAPQ:
798 defaultMapq = parseInt(0, "--mapq must be positive");
799 break;
800 case ARG_MAXBTS: {
801 maxBts = parseInt(0, "--maxbts must be positive");
802 maxBtsBetter = maxBts;
803 break;
804 }
805 case ARG_STRAND_FIX: strandFix = true; break;
806 case ARG_PARTITION: partitionSz = parse<int>(optarg); break;
807 case ARG_READS_PER_BATCH: {
808 if(optarg == NULL || parse<int>(optarg) < 1) {
809 cerr << "--reads-per-batch arg must be at least 1" << endl;
810 printUsage(cerr);
811 throw 1;
812 }
813 // TODO: should output batch size be controlled separately?
814 readsPerBatch = outBatchSz = parse<int>(optarg);
815 break;
816 }
817 case ARG_ORIG:
818 if(optarg == NULL || strlen(optarg) == 0) {
819 cerr << "--orig arg must be followed by a string" << endl;
820 printUsage(cerr);
821 throw 1;
822 }
823 origString = optarg;
824 break;
825 case -1: break; /* Done with options. */
826 case 0:
827 if (long_options[option_index].flag != 0)
828 break;
829 default:
830 printUsage(cerr);
831 throw 1;
832 }
833 } while(next_option != -1);
834
835 if (nthreads == 1 && !thread_stealing) {
836 reorder = false;
837 }
838 if (reorder == true && outType != OUTPUT_SAM) {
839 cerr << "Bowtie will reorder its output only when outputting SAM." << endl
840 << "Please specify the `-S` parameter if you intend on using this option." << endl;
841 throw 1;
842 }
843 //bool paired = mates1.size() > 0 || mates2.size() > 0 || mates12.size() > 0;
844 if(rangeMode) {
845 // Tell the Ebwt loader to ignore the suffix-array portion of
846 // the index. We don't need it because the user isn't asking
847 // for bowtie to report reference positions (just matrix
848 // ranges).
849 offRate = 32;
850 }
851 if(!maqLike && mismatches == 3) {
852 // Much faster than normal 3-mismatch mode
853 stateful = true;
854 }
855 if(mates1.size() != mates2.size()) {
856 cerr << "Error: " << mates1.size() << " mate files/sequences were specified with -1, but " << mates2.size() << endl
857 << "mate files/sequences were specified with -2. The same number of mate files/" << endl
858 << "sequences must be specified with -1 and -2." << endl;
859 throw 1;
860 }
861 // Check for duplicate mate input files
862 if(format != CMDLINE) {
863 for(size_t i = 0; i < mates1.size(); i++) {
864 for(size_t j = 0; j < mates2.size(); j++) {
865 if(mates1[i] == mates2[j] && !quiet) {
866 cerr << "Warning: Same mate file \"" << mates1[i] << "\" appears as argument to both -1 and -2" << endl;
867 }
868 }
869 }
870 }
871 if(tryHard) {
872 // Increase backtracking limit to huge number
873 maxBts = maxBtsBetter = INT_MAX;
874 // Increase number of paired-end scan attempts to huge number
875 mixedAttemptLim = UINT_MAX;
876 }
877 if(!stateful && sampleMax) {
878 if(!quiet) {
879 cerr << "Warning: -M was specified w/o --best; automatically enabling --best" << endl;
880 }
881 stateful = true;
882 }
883 if(strata && !stateful) {
884 cerr << "--strata must be combined with --best" << endl;
885 throw 1;
886 }
887 if(strata && !allHits && khits == 1 && mhits == 0xffffffff) {
888 cerr << "--strata has no effect unless combined with -m, -a, or -k N where N > 1" << endl;
889 throw 1;
890 }
891 // If both -s and -u are used, we need to adjust qUpto accordingly
892 // since it uses patid to know if we've reached the -u limit (and
893 // patids are all shifted up by skipReads characters)
894 if(qUpto + skipReads > qUpto) {
895 qUpto += skipReads;
896 }
897 if(useShmem && useMm && !quiet) {
898 cerr << "Warning: --shmem overrides --mm..." << endl;
899 useMm = false;
900 }
901 if(!mateFwSet) {
902 // Set nucleotide space default (--fr)
903 mate1fw = true;
904 mate2fw = false;
905 }
906 if(outType != OUTPUT_FULL && suppressOuts.count() > 0 && !quiet) {
907 cerr << "Warning: Ignoring --suppress because output type is not default." << endl;
908 cerr << " --suppress is only available for the default output type." << endl;
909 suppressOuts.clear();
910 }
911 thread_stealing = thread_ceiling > nthreads;
912 #ifdef _WIN32
913 thread_stealing = false;
914 #endif
915 if(thread_stealing && thread_stealing_dir.empty()) {
916 cerr << "When --thread-ceiling is specified, must also specify --thread-piddir" << endl;
917 throw 1;
918 }
919 }
920
921 static const char *argv0 = NULL;
922
923 #define FINISH_READ(p) \
924 /* Don't do finishRead if the read isn't legit or if the read was skipped by the doneMask */ \
925 if(get_read_ret.first) { \
926 sink->finishRead(*p, true, !skipped); \
927 get_read_ret.first = false; \
928 } \
929 skipped = false;
930
931 /// Macro for getting the next read, possibly aborting depending on
932 /// whether the result is empty or the patid exceeds the limit, and
933 /// marshaling the read into convenient variables.
934 #define GET_READ(p) \
935 if(get_read_ret.second) break; \
936 get_read_ret = p->nextReadPair(); \
937 if(p->rdid() >= qUpto) { \
938 get_read_ret = make_pair(false, true); \
939 } \
940 if(!get_read_ret.first) { \
941 if(get_read_ret.second) { \
942 break; \
943 } \
944 continue; \
945 } \
946 BTDnaString& patFw = p->bufa().patFw; \
947 patFw.length(); \
948 BTDnaString& patRc = p->bufa().patRc; \
949 patRc.length(); \
950 BTString& qual = p->bufa().qual; \
951 qual.length(); \
952 BTString& qualRev = p->bufa().qualRev; \
953 qualRev.length(); \
954 BTDnaString& patFwRev = p->bufa().patFwRev; \
955 patFwRev.length(); \
956 BTDnaString& patRcRev = p->bufa().patRcRev; \
957 patRcRev.length(); \
958 BTString& name = p->bufa().name; \
959 name.length(); \
960 uint32_t patid = (uint32_t)p->rdid(); \
961 params.setPatId(patid);
962
963 #define WORKER_EXIT() \
964 patsrcFact->destroy(patsrc); \
965 delete patsrcFact; \
966 sinkFact->destroy(sink); \
967 delete sinkFact; \
968 return;
969
970 #ifdef CHUD_PROFILING
971 #define CHUD_START() chudStartRemotePerfMonitor("Bowtie");
972 #define CHUD_STOP() chudStopRemotePerfMonitor();
973 #else
974 #define CHUD_START()
975 #define CHUD_STOP()
976 #endif
977
978 /// Create a PatternSourcePerThread for the current thread according
979 /// to the global params and return a pointer to it
980 static PatternSourcePerThreadFactory*
createPatsrcFactory(PatternComposer & _patsrc,int tid,uint32_t max_buf)981 createPatsrcFactory(PatternComposer& _patsrc, int tid, uint32_t max_buf) {
982 PatternSourcePerThreadFactory *patsrcFact;
983 patsrcFact = new PatternSourcePerThreadFactory(_patsrc, max_buf, skipReads, seed);
984 assert(patsrcFact != NULL);
985 return patsrcFact;
986 }
987
988 /**
989 * Allocate a HitSinkPerThreadFactory on the heap according to the
990 * global params and return a pointer to it.
991 */
992 static HitSinkPerThreadFactory*
createSinkFactory(HitSink & _sink,size_t threadId)993 createSinkFactory(HitSink& _sink, size_t threadId) {
994 HitSinkPerThreadFactory *sink = NULL;
995 if(!strata) {
996 // Unstratified
997 if(!allHits) {
998 // First N good; "good" inherently ignores strata
999 sink = new NGoodHitSinkPerThreadFactory(_sink, khits, mhits, defaultMapq, threadId);
1000 } else {
1001 // All hits, spanning strata
1002 sink = new AllHitSinkPerThreadFactory(_sink, mhits, defaultMapq, threadId);
1003 }
1004 } else {
1005 // Stratified
1006 assert(stateful);
1007 if(!allHits) {
1008 assert(stateful);
1009 // Buffer best hits, assuming they're arriving in best-
1010 // to-worst order
1011 sink = new NBestFirstStratHitSinkPerThreadFactory(_sink, khits, mhits, defaultMapq, threadId);
1012 } else {
1013 assert(stateful);
1014 // Buffer best hits, assuming they're arriving in best-
1015 // to-worst order
1016 sink = new NBestFirstStratHitSinkPerThreadFactory(_sink, 0xffffffff/2, mhits, defaultMapq, threadId);
1017 }
1018 }
1019 assert(sink != NULL);
1020 return sink;
1021 }
1022
increment_thread_counter()1023 void increment_thread_counter() {
1024 #if (__cplusplus >= 201103)
1025 thread_counter.fetch_add(1);
1026 #else
1027 ThreadSafe ts(&thread_counter_mutex);
1028 thread_counter++;
1029 #endif
1030 }
1031
decrement_thread_counter()1032 void decrement_thread_counter() {
1033 #if (__cplusplus >= 201103)
1034 thread_counter.fetch_sub(1);
1035 #else
1036 ThreadSafe ts(&thread_counter_mutex);
1037 thread_counter--;
1038 #endif
1039 }
1040
1041 #ifndef _WIN32
del_pid(const char * dirname,int pid)1042 void del_pid(const char* dirname,int pid) {
1043 struct stat finfo;
1044 char* fname = (char*) calloc(FNAME_SIZE,sizeof(char));
1045 sprintf(fname,"%s/%d",dirname,pid);
1046 if(stat( fname, &finfo) != 0) {
1047 free(fname);
1048 return;
1049 }
1050 unlink(fname);
1051 free(fname);
1052 }
1053
1054 //from http://stackoverflow.com/questions/18100097/portable-way-to-check-if-directory-exists-windows-linux-c
write_pid(const char * dirname,int pid)1055 static void write_pid(const char* dirname,int pid) {
1056 struct stat dinfo;
1057 if(stat(dirname, &dinfo) != 0) {
1058 mkdir(dirname,0755);
1059 }
1060 char* fname = (char*) calloc(FNAME_SIZE,sizeof(char));
1061 sprintf(fname,"%s/%d",dirname,pid);
1062 FILE* f = fopen(fname,"w");
1063 fclose(f);
1064 free(fname);
1065 }
1066
1067 //from http://stackoverflow.com/questions/612097/how-can-i-get-the-list-of-files-in-a-directory-using-c-or-c
read_dir(const char * dirname,int * num_pids)1068 static int read_dir(const char* dirname,int* num_pids) {
1069 DIR *dir;
1070 struct dirent *ent;
1071 char* fname = (char*) calloc(FNAME_SIZE,sizeof(char));
1072 int lowest_pid = -1;
1073 if ((dir = opendir (dirname)) != NULL) {
1074 while ((ent = readdir (dir)) != NULL) {
1075 if(ent->d_name[0] == '.')
1076 continue;
1077 int pid = atoi(ent->d_name);
1078 sprintf(fname,"/proc/%s", ent->d_name);
1079 if(kill(pid, 0) != 0) {
1080 //deleting pids can lead to race conditions if
1081 //2 or more BT2 processes both try to delete
1082 //so just skip instead
1083 //del_pid(dirname,pid);
1084 continue;
1085 }
1086 (*num_pids)++;
1087 if(pid < lowest_pid || lowest_pid == -1)
1088 lowest_pid = pid;
1089 }
1090 closedir (dir);
1091 } else {
1092 perror (""); // could not open directory
1093 }
1094 free(fname);
1095 return lowest_pid;
1096 }
1097
steal_thread(int pid,int orig_nthreads)1098 static bool steal_thread(int pid, int orig_nthreads) {
1099 int ncpu = thread_ceiling;
1100 if(thread_ceiling <= nthreads) {
1101 return false;
1102 }
1103 int num_pids = 0;
1104 int lowest_pid = read_dir(thread_stealing_dir.c_str(), &num_pids);
1105 if(lowest_pid != pid) {
1106 return false;
1107 }
1108 int in_use = ((num_pids-1) * orig_nthreads) + nthreads; //in_use is now baseline + ours
1109 float spare = ncpu - in_use;
1110 int spare_r = floor(spare);
1111 float r = rand() % 100/100.0;
1112 if(r <= (spare - spare_r)) {
1113 spare_r = ceil(spare);
1114 }
1115 return spare_r > 0;
1116 }
1117 #endif
1118
1119 /**
1120 * Search through a single (forward) Ebwt index for exact end-to-end
1121 * hits. Assumes that index is already loaded into memory.
1122 */
1123 static PatternComposer* exactSearch_patsrc;
1124 static HitSink* exactSearch_sink;
1125 static Ebwt* exactSearch_ebwt;
1126 static EList<BTRefString>* exactSearch_os;
1127 static BitPairReference* exactSearch_refs;
1128 #if (__cplusplus >= 201103L)
1129 //void exactSearchWorker::operator()() const {
exactSearchWorker(void * vp)1130 static void exactSearchWorker(void *vp) {
1131 thread_tracking_pair *p = (thread_tracking_pair*) vp;
1132 int tid = p->tid;
1133 #else
1134 static void exactSearchWorker(void *vp) {
1135 int tid = *((int*)vp);
1136 #endif
1137 if(thread_stealing) {
1138 increment_thread_counter();
1139 }
1140 PatternComposer& _patsrc = *exactSearch_patsrc;
1141 HitSink& _sink = *exactSearch_sink;
1142 Ebwt& ebwt = *exactSearch_ebwt;
1143 EList<BTRefString >& os = *exactSearch_os;
1144
1145 // Per-thread initialization
1146 PatternSourcePerThreadFactory *patsrcFact = createPatsrcFactory(_patsrc, tid, readsPerBatch);
1147 PatternSourcePerThread *patsrc = patsrcFact->create();
1148 HitSinkPerThreadFactory* sinkFact = createSinkFactory(_sink, tid);
1149 HitSinkPerThread* sink = sinkFact->create();
1150 EbwtSearchParams params(
1151 *sink, // HitSink
1152 os, // reference sequences
1153 true, // read is forward
1154 true); // index is forward
1155 GreedyDFSRangeSource bt(
1156 &ebwt, params,
1157 0xffffffff, // qualThresh
1158 0xffffffff, // max backtracks (no max)
1159 0, // reportPartials (don't)
1160 true, // reportExacts
1161 rangeMode, // reportRanges
1162 NULL, // seedlings
1163 NULL, // mutations
1164 verbose, // verbose
1165 &os,
1166 false); // considerQuals
1167 pair<bool, bool> get_read_ret = make_pair(false, false);
1168 bool skipped = false;
1169 #ifdef PER_THREAD_TIMING
1170 uint64_t ncpu_changeovers = 0;
1171 uint64_t nnuma_changeovers = 0;
1172
1173 int current_cpu = 0, current_node = 0;
1174 get_cpu_and_node(current_cpu, current_node);
1175
1176 std::stringstream ss;
1177 std::string msg;
1178 ss << "thread: " << tid << " time: ";
1179 msg = ss.str();
1180 {
1181 Timer timer(std::cout, msg.c_str());
1182 #endif
1183 while(true) {
1184 #ifdef PER_THREAD_TIMING
1185 int cpu = 0, node = 0;
1186 get_cpu_and_node(cpu, node);
1187 if(cpu != current_cpu) {
1188 ncpu_changeovers++;
1189 current_cpu = cpu;
1190 }
1191 if(node != current_node) {
1192 nnuma_changeovers++;
1193 current_node = node;
1194 }
1195 #endif
1196 FINISH_READ(patsrc);
1197 GET_READ(patsrc);
1198 #include "search_exact.c"
1199 }
1200 FINISH_READ(patsrc);
1201 if(thread_stealing) {
1202 decrement_thread_counter();
1203 }
1204 #ifdef PER_THREAD_TIMING
1205 ss.str("");
1206 ss.clear();
1207 ss << "thread: " << tid << " cpu_changeovers: " << ncpu_changeovers << std::endl
1208 << "thread: " << tid << " node_changeovers: " << nnuma_changeovers << std::endl;
1209 std::cout << ss.str();
1210 }
1211 #endif
1212 #if (__cplusplus >= 201103L)
1213 p->done->fetch_add(1);
1214 #endif
1215 WORKER_EXIT();
1216 }
1217
1218 /**
1219 * A statefulness-aware worker driver. Uses UnpairedExactAlignerV1.
1220 */
1221 #if (__cplusplus >= 201103L)
1222 //void exactSearchWorkerStateful::operator()() const {
1223 static void exactSearchWorkerStateful(void *vp) {
1224 thread_tracking_pair *p = (thread_tracking_pair*) vp;
1225 int tid = p->tid;
1226 #else
1227 static void exactSearchWorkerStateful(void *vp) {
1228 int tid = *((int*)vp);
1229 #endif
1230 if(thread_stealing) {
1231 increment_thread_counter();
1232 }
1233 PatternComposer& _patsrc = *exactSearch_patsrc;
1234 HitSink& _sink = *exactSearch_sink;
1235 Ebwt& ebwt = *exactSearch_ebwt;
1236 EList<BTRefString >& os = *exactSearch_os;
1237 BitPairReference* refs = exactSearch_refs;
1238
1239 // Global initialization
1240 PatternSourcePerThreadFactory* patsrcFact = createPatsrcFactory(_patsrc, tid, readsPerBatch);
1241 HitSinkPerThreadFactory* sinkFact = createSinkFactory(_sink, tid);
1242
1243 ChunkPool *pool = new ChunkPool(chunkSz * 1024, chunkPoolMegabytes * 1024 * 1024, chunkVerbose);
1244 UnpairedExactAlignerV1Factory alSEfact(
1245 ebwt,
1246 !nofw,
1247 !norc,
1248 _sink,
1249 *sinkFact,
1250 NULL, //&cacheFw,
1251 NULL, //&cacheBw,
1252 cacheLimit,
1253 pool,
1254 refs,
1255 os,
1256 !noMaqRound,
1257 !better,
1258 strandFix,
1259 rangeMode,
1260 verbose,
1261 quiet,
1262 seed);
1263 PairedExactAlignerV1Factory alPEfact(
1264 ebwt,
1265 !nofw,
1266 !norc,
1267 useV1,
1268 _sink,
1269 *sinkFact,
1270 mate1fw,
1271 mate2fw,
1272 minInsert,
1273 maxInsert,
1274 dontReconcileMates,
1275 mhits, // for symCeiling
1276 mixedThresh,
1277 mixedAttemptLim,
1278 NULL, //&cacheFw,
1279 NULL, //&cacheBw,
1280 cacheLimit,
1281 pool,
1282 refs, os,
1283 reportSe,
1284 !noMaqRound,
1285 strandFix,
1286 !better,
1287 rangeMode,
1288 verbose,
1289 quiet,
1290 seed);
1291 {
1292 MixedMultiAligner multi(
1293 prefetchWidth,
1294 qUpto,
1295 alSEfact,
1296 alPEfact,
1297 *patsrcFact);
1298 // Run that mother
1299 multi.run(false, tid);
1300 // MultiAligner must be destroyed before patsrcFact
1301 }
1302
1303 if(thread_stealing) {
1304 decrement_thread_counter();
1305 }
1306
1307 delete patsrcFact;
1308 delete sinkFact;
1309 delete pool;
1310 #if (__cplusplus >= 201103L)
1311 p->done->fetch_add(1);
1312 #endif
1313 return;
1314 }
1315
1316 #define SET_A_FW(bt, p, params) \
1317 bt.setQuery(&p->bufa().patFw, &p->bufa().qual, &p->bufa().name); \
1318 params.setFw(true);
1319 #define SET_A_RC(bt, p, params) \
1320 bt.setQuery(&p->bufa().patRc, &p->bufa().qualRev, &p->bufa().name); \
1321 params.setFw(false);
1322 #define SET_B_FW(bt, p, params) \
1323 bt.setQuery(&p->bufb().patFw, &p->bufb().qual, &p->bufb().name); \
1324 params.setFw(true);
1325 #define SET_B_RC(bt, p, params) \
1326 bt.setQuery(&p->bufb().patRc, &p->bufb().qualRev, &p->bufb().name); \
1327 params.setFw(false);
1328
1329 /**
1330 * Search through a single (forward) Ebwt index for exact end-to-end
1331 * hits. Assumes that index is already loaded into memory.
1332 */
1333 static void exactSearch(PatternComposer& _patsrc,
1334 HitSink& _sink,
1335 Ebwt& ebwt,
1336 EList<BTRefString >& os)
1337 {
1338 exactSearch_patsrc = &_patsrc;
1339 exactSearch_sink = &_sink;
1340 exactSearch_ebwt = &ebwt;
1341 exactSearch_os = &os;
1342
1343 assert(!ebwt.isInMemory());
1344 {
1345 // Load the rest of (vast majority of) the backward Ebwt into
1346 // memory
1347 Timer _t(cerr, "Time loading forward index: ", timing);
1348 ebwt.loadIntoMemory(-1, !noRefNames, startVerbose);
1349 }
1350
1351 BitPairReference *refs = NULL;
1352 bool pair = mates1.size() > 0 || mates12.size() > 0;
1353 if((pair && mixedThresh < 0xffffffff)) {
1354 Timer _t(cerr, "Time loading reference: ", timing);
1355 refs = new BitPairReference(adjustedEbwtFileBase, sanityCheck, NULL, &os, false, true, useMm, useShmem, mmSweep, verbose, startVerbose);
1356 if(!refs->loaded()) throw 1;
1357 }
1358 exactSearch_refs = refs;
1359 int tids[max(nthreads, thread_ceiling)];
1360 #if (__cplusplus >= 201103L)
1361 EList<std::thread*> threads;
1362 EList<thread_tracking_pair> tps;
1363 tps.resizeExact(max(nthreads, thread_ceiling));
1364 #else
1365 EList<tthread::thread*> threads;
1366 #endif
1367
1368 #if (__cplusplus >= 201103L)
1369 std::atomic<int> all_threads_done;
1370 all_threads_done = 0;
1371 #endif
1372 CHUD_START();
1373 {
1374 Timer _t(cerr, "Time for 0-mismatch search: ", timing);
1375
1376 #ifndef _WIN32
1377 int pid = 0;
1378 if(thread_stealing) {
1379 pid = getpid();
1380 write_pid(thread_stealing_dir.c_str(), pid);
1381 thread_counter = 0;
1382 }
1383 #endif
1384
1385 for(int i = 0; i < nthreads; i++) {
1386 tids[i] = i;
1387 #if (__cplusplus >= 201103L)
1388 tps[i].tid = tids[i];
1389 tps[i].done = &all_threads_done;
1390 if (i == nthreads - 1) {
1391 if(stateful) {
1392 exactSearchWorkerStateful((void*)&tps[i]);
1393 } else {
1394 exactSearchWorker((void*)&tps[i]);
1395 }
1396 }
1397 else {
1398 if(stateful) {
1399 threads.push_back(new std::thread(exactSearchWorkerStateful, (void*)&tps[i]));
1400 } else {
1401 threads.push_back(new std::thread(exactSearchWorker, (void*)&tps[i]));
1402 }
1403 threads[i]->detach();
1404 SLEEP(10);
1405 }
1406 #else
1407 if (i == nthreads - 1) {
1408 if(stateful) {
1409 exactSearchWorkerStateful((void*)(tids + i));
1410 } else {
1411 exactSearchWorker((void*)(tids + i));
1412 }
1413 }
1414 else {
1415 if(stateful) {
1416 threads.push_back(new tthread::thread(exactSearchWorkerStateful, (void *)(tids + i)));
1417 } else {
1418 threads.push_back(new tthread::thread(exactSearchWorker, (void *)(tids + i)));
1419 }
1420 }
1421 #endif
1422 }
1423
1424 #ifndef _WIN32
1425 if(thread_stealing) {
1426 int orig_threads = nthreads, steal_ctr = 1;
1427 for(int j = 0; j < 10; j++) {
1428 sleep(1);
1429 }
1430 while(thread_counter > 0) {
1431 if(steal_thread(pid, orig_threads)) {
1432 nthreads++;
1433 tids[nthreads-1] = nthreads;
1434 #if (__cplusplus >= 201103L)
1435 tps[nthreads-1].tid = nthreads - 1;
1436 tps[nthreads-1].done = &all_threads_done;
1437 if(stateful) {
1438 threads.push_back(new std::thread(exactSearchWorkerStateful, (void*)&tps[nthreads-1]));
1439 } else {
1440 threads.push_back(new std::thread(exactSearchWorker, (void*)&tps[nthreads-1]));
1441 }
1442 threads[nthreads-1]->detach();
1443 SLEEP(10);
1444 #else
1445 if(stateful) {
1446 threads.push_back(new tthread::thread(exactSearchWorkerStateful, (void *)(tids + nthreads - 1)));
1447 } else {
1448 threads.push_back(new tthread::thread(exactSearchWorker, (void *)(tids + nthreads - 1)));
1449 }
1450 #endif
1451 cerr << "pid " << pid << " started new worker # " << nthreads << endl;
1452 }
1453 steal_ctr++;
1454 for(int j = 0; j < 10; j++) {
1455 sleep(1);
1456 }
1457 }
1458 }
1459 #endif
1460
1461 #if (__cplusplus >= 201103L)
1462 while(all_threads_done < nthreads) {
1463 SLEEP(10);
1464 }
1465 #else
1466 for (int i = 0; i < nthreads - 1; i++) {
1467 threads[i]->join();
1468 }
1469 #endif
1470
1471 #ifndef _WIN32
1472 if(thread_stealing) {
1473 del_pid(thread_stealing_dir.c_str(), pid);
1474 }
1475 #endif
1476 }
1477 if(refs != NULL) delete refs;
1478
1479 for (int i = 0; i < nthreads - 1; i++) {
1480 if (threads[i] != NULL) {
1481 delete threads[i];
1482 }
1483 }
1484 }
1485
1486 /**
1487 * Search through a pair of Ebwt indexes, one for the forward direction
1488 * and one for the backward direction, for exact end-to-end hits and 1-
1489 * mismatch end-to-end hits. In my experience, this is slightly faster
1490 * than Maq (default) mode with the -n 1 option.
1491 *
1492 * Forward Ebwt (ebwtFw) is already loaded into memory and backward
1493 * Ebwt (ebwtBw) is not loaded into memory.
1494 */
1495 static PatternComposer* mismatchSearch_patsrc;
1496 static HitSink* mismatchSearch_sink;
1497 static Ebwt* mismatchSearch_ebwtFw;
1498 static Ebwt* mismatchSearch_ebwtBw;
1499 static EList<BTRefString >* mismatchSearch_os;
1500 static SyncBitset* mismatchSearch_doneMask;
1501 static SyncBitset* mismatchSearch_hitMask;
1502 static BitPairReference* mismatchSearch_refs;
1503
1504 /**
1505 * A statefulness-aware worker driver. Uses Unpaired/Paired1mmAlignerV1.
1506 */
1507 #if (__cplusplus >= 201103L)
1508 //void mismatchSearchWorkerFullStateful::operator()() const {
1509 static void mismatchSearchWorkerFullStateful(void *vp) {
1510 thread_tracking_pair *p = (thread_tracking_pair*) vp;
1511 int tid = p->tid;
1512 #else
1513 static void mismatchSearchWorkerFullStateful(void *vp) {
1514 int tid = *((int*)vp);
1515 #endif
1516 if(thread_stealing) {
1517 increment_thread_counter();
1518 }
1519 PatternComposer& _patsrc = *mismatchSearch_patsrc;
1520 HitSink& _sink = *mismatchSearch_sink;
1521 Ebwt& ebwtFw = *mismatchSearch_ebwtFw;
1522 Ebwt& ebwtBw = *mismatchSearch_ebwtBw;
1523 EList<BTRefString >& os = *mismatchSearch_os;
1524 BitPairReference* refs = mismatchSearch_refs;
1525
1526 // Global initialization
1527 PatternSourcePerThreadFactory* patsrcFact = createPatsrcFactory(_patsrc, tid, readsPerBatch);
1528 HitSinkPerThreadFactory* sinkFact = createSinkFactory(_sink, tid);
1529 ChunkPool *pool = new ChunkPool(chunkSz * 1024, chunkPoolMegabytes * 1024 * 1024, chunkVerbose);
1530
1531 Unpaired1mmAlignerV1Factory alSEfact(
1532 ebwtFw,
1533 &ebwtBw,
1534 !nofw,
1535 !norc,
1536 _sink,
1537 *sinkFact,
1538 NULL, //&cacheFw,
1539 NULL, //&cacheBw,
1540 cacheLimit,
1541 pool,
1542 refs,
1543 os,
1544 !noMaqRound,
1545 !better,
1546 strandFix,
1547 rangeMode,
1548 verbose,
1549 quiet,
1550 seed);
1551 Paired1mmAlignerV1Factory alPEfact(
1552 ebwtFw,
1553 &ebwtBw,
1554 !nofw,
1555 !norc,
1556 useV1,
1557 _sink,
1558 *sinkFact,
1559 mate1fw,
1560 mate2fw,
1561 minInsert,
1562 maxInsert,
1563 dontReconcileMates,
1564 mhits, // for symCeiling
1565 mixedThresh,
1566 mixedAttemptLim,
1567 NULL, //&cacheFw,
1568 NULL, //&cacheBw,
1569 cacheLimit,
1570 pool,
1571 refs, os,
1572 reportSe,
1573 !noMaqRound,
1574 !better,
1575 strandFix,
1576 rangeMode,
1577 verbose,
1578 quiet,
1579 seed);
1580 {
1581 MixedMultiAligner multi(
1582 prefetchWidth,
1583 qUpto,
1584 alSEfact,
1585 alPEfact,
1586 *patsrcFact);
1587 // Run that mother
1588 multi.run(false, tid);
1589 // MultiAligner must be destroyed before patsrcFact
1590 }
1591
1592 if(thread_stealing) {
1593 decrement_thread_counter();
1594 }
1595 #if (__cplusplus >= 201103L)
1596 p->done->fetch_add(1);
1597 #endif
1598
1599 delete patsrcFact;
1600 delete sinkFact;
1601 delete pool;
1602 return;
1603 }
1604 #if (__cplusplus >= 201103L)
1605 //void mismatchSearchWorkerFull::operator()() const {
1606 static void mismatchSearchWorkerFull(void *vp){
1607 thread_tracking_pair *p = (thread_tracking_pair*) vp;
1608 int tid = p->tid;
1609 #else
1610 static void mismatchSearchWorkerFull(void *vp){
1611 int tid = *((int*)vp);
1612 #endif
1613 if(thread_stealing) {
1614 increment_thread_counter();
1615 }
1616 PatternComposer& _patsrc = *mismatchSearch_patsrc;
1617 HitSink& _sink = *mismatchSearch_sink;
1618 Ebwt& ebwtFw = *mismatchSearch_ebwtFw;
1619 Ebwt& ebwtBw = *mismatchSearch_ebwtBw;
1620 EList<BTRefString >& os = *mismatchSearch_os;
1621
1622 // Per-thread initialization
1623 PatternSourcePerThreadFactory* patsrcFact = createPatsrcFactory(_patsrc, tid, readsPerBatch);
1624 PatternSourcePerThread* patsrc = patsrcFact->create();
1625 HitSinkPerThreadFactory* sinkFact = createSinkFactory(_sink, tid);
1626 HitSinkPerThread* sink = sinkFact->create();
1627 EbwtSearchParams params(
1628 *sink, // HitSinkPerThread
1629 os, // reference sequences
1630 true, // read is forward
1631 false); // index is mirror index
1632 GreedyDFSRangeSource bt(
1633 &ebwtFw, params,
1634 0xffffffff, // qualThresh
1635 0xffffffff, // max backtracks (no max)
1636 0, // reportPartials (don't)
1637 true, // reportExacts
1638 rangeMode, // reportRanges
1639 NULL, // seedlings
1640 NULL, // mutations
1641 verbose, // verbose
1642 &os,
1643 false); // considerQuals
1644 bool skipped = false;
1645 pair<bool, bool> get_read_ret = make_pair(false, false);
1646 #ifdef PER_THREAD_TIMING
1647 uint64_t ncpu_changeovers = 0;
1648 uint64_t nnuma_changeovers = 0;
1649
1650 int current_cpu = 0, current_node = 0;
1651 get_cpu_and_node(current_cpu, current_node);
1652
1653 std::stringstream ss;
1654 std::string msg;
1655 ss << "thread: " << tid << " time: ";
1656 msg = ss.str();
1657 {
1658 Timer timer(std::cout, msg.c_str());
1659 #endif
1660 while(true) {
1661 #ifdef PER_THREAD_TIMING
1662 int cpu = 0, node = 0;
1663 get_cpu_and_node(cpu, node);
1664 if(cpu != current_cpu) {
1665 ncpu_changeovers++;
1666 current_cpu = cpu;
1667 }
1668 if(node != current_node) {
1669 nnuma_changeovers++;
1670 current_node = node;
1671 }
1672 #endif
1673 FINISH_READ(patsrc);
1674 GET_READ(patsrc);
1675 uint32_t plen = (uint32_t)patFw.length();
1676 uint32_t s = plen;
1677 uint32_t s3 = s >> 1; // length of 3' half of seed
1678 uint32_t s5 = (s >> 1) + (s & 1); // length of 5' half of seed
1679 #define DONEMASK_SET(p)
1680 #include "search_1mm_phase1.c"
1681 #include "search_1mm_phase2.c"
1682 #undef DONEMASK_SET
1683 } // End read loop
1684 FINISH_READ(patsrc);
1685 #ifdef PER_THREAD_TIMING
1686 ss.str("");
1687 ss.clear();
1688 ss << "thread: " << tid << " cpu_changeovers: " << ncpu_changeovers << std::endl
1689 << "thread: " << tid << " node_changeovers: " << nnuma_changeovers << std::endl;
1690 std::cout << ss.str();
1691 }
1692 #endif
1693 #if (__cplusplus >= 201103L)
1694 p->done->fetch_add(1);
1695 #endif
1696 WORKER_EXIT();
1697 }
1698
1699 /**
1700 * Search through a single (forward) Ebwt index for exact end-to-end
1701 * hits. Assumes that index is already loaded into memory.
1702 */
1703 static void mismatchSearchFull(PatternComposer& _patsrc,
1704 HitSink& _sink,
1705 Ebwt& ebwtFw,
1706 Ebwt& ebwtBw,
1707 EList<BTRefString >& os)
1708 {
1709 mismatchSearch_patsrc = &_patsrc;
1710 mismatchSearch_sink = &_sink;
1711 mismatchSearch_ebwtFw = &ebwtFw;
1712 mismatchSearch_ebwtBw = &ebwtBw;
1713 mismatchSearch_doneMask = NULL;
1714 mismatchSearch_hitMask = NULL;
1715 mismatchSearch_os = &os;
1716
1717 assert(!ebwtFw.isInMemory());
1718 assert(!ebwtBw.isInMemory());
1719 {
1720 // Load the other half of the index into memory
1721 Timer _t(cerr, "Time loading forward index: ", timing);
1722 ebwtFw.loadIntoMemory(-1, !noRefNames, startVerbose);
1723 }
1724 {
1725 // Load the other half of the index into memory
1726 Timer _t(cerr, "Time loading mirror index: ", timing);
1727 ebwtBw.loadIntoMemory(-1, !noRefNames, startVerbose);
1728 }
1729 // Create range caches, which are shared among all aligners
1730 BitPairReference *refs = NULL;
1731 bool pair = mates1.size() > 0 || mates12.size() > 0;
1732 if(pair && mixedThresh < 0xffffffff) {
1733 Timer _t(cerr, "Time loading reference: ", timing);
1734 refs = new BitPairReference(adjustedEbwtFileBase, sanityCheck, NULL, &os, false, true, useMm, useShmem, mmSweep, verbose, startVerbose);
1735 if(!refs->loaded()) throw 1;
1736 }
1737 mismatchSearch_refs = refs;
1738
1739 int tids[max(nthreads, thread_ceiling)];
1740 #if (__cplusplus >= 201103L)
1741 EList<std::thread*> threads;
1742 EList<thread_tracking_pair> tps;
1743 tps.resizeExact(max(nthreads, thread_ceiling));
1744 #else
1745 EList<tthread::thread*> threads;
1746 #endif
1747
1748 #if (__cplusplus >= 201103L)
1749 std::atomic<int> all_threads_done;
1750 all_threads_done = 0;
1751 #endif
1752
1753 CHUD_START();
1754 {
1755 Timer _t(cerr, "Time for 1-mismatch full-index search: ", timing);
1756
1757 #ifndef _WIN32
1758 int pid = 0;
1759 if(thread_stealing) {
1760 pid = getpid();
1761 write_pid(thread_stealing_dir.c_str(), pid);
1762 thread_counter = 0;
1763 }
1764 #endif
1765
1766 for(int i = 0; i < nthreads; i++) {
1767 tids[i] = i;
1768 #if (__cplusplus >= 201103L)
1769 tps[i].tid = tids[i];
1770 tps[i].done = &all_threads_done;
1771 if (i == nthreads - 1) {
1772 if(stateful) {
1773 mismatchSearchWorkerFullStateful((void*)&tps[i]);
1774 } else {
1775 mismatchSearchWorkerFull((void*)&tps[i]);
1776 }
1777 }
1778 else {
1779 if(stateful) {
1780 threads.push_back(new std::thread(mismatchSearchWorkerFullStateful, (void*)&tps[i]));
1781 } else {
1782 threads.push_back(new std::thread(mismatchSearchWorkerFull, (void*)&tps[i]));
1783 }
1784 threads[i]->detach();
1785 SLEEP(10);
1786 }
1787 #else
1788 if (i == nthreads - 1) {
1789 if(stateful) {
1790 mismatchSearchWorkerFullStateful((void*)(tids + i));
1791 } else {
1792 mismatchSearchWorkerFull((void*)(tids + i));
1793 }
1794 }
1795 else {
1796 if(stateful) {
1797 threads.push_back(new tthread::thread(mismatchSearchWorkerFullStateful, (void *)(tids + i)));
1798 } else {
1799 threads.push_back(new tthread::thread(mismatchSearchWorkerFull, (void *)(tids + i)));
1800 }
1801 }
1802 #endif
1803 }
1804
1805 #ifndef _WIN32
1806 if(thread_stealing) {
1807 int orig_threads = nthreads, steal_ctr = 1;
1808 for(int j = 0; j < 10; j++) {
1809 sleep(1);
1810 }
1811 while(thread_counter > 0) {
1812 if(steal_thread(pid, orig_threads)) {
1813 nthreads++;
1814 tids[nthreads-1] = nthreads;
1815 #if (__cplusplus >= 201103L)
1816 tps[nthreads-1].tid = nthreads - 1;
1817 tps[nthreads-1].done = &all_threads_done;
1818 if(stateful) {
1819 threads.push_back(new std::thread(mismatchSearchWorkerFullStateful, (void*)&tps[nthreads-1]));
1820 } else {
1821 threads.push_back(new std::thread(mismatchSearchWorkerFull, (void*)&tps[nthreads-1]));
1822 }
1823 threads[nthreads - 1]->detach();
1824 SLEEP(10);
1825 #else
1826 if(stateful) {
1827 threads.push_back(new tthread::thread(mismatchSearchWorkerFullStateful, (void *)(tids + nthreads - 1)));
1828 } else {
1829 threads.push_back(new tthread::thread(mismatchSearchWorkerFull, (void *)(tids + nthreads - 1)));
1830 }
1831 #endif
1832 cerr << "pid " << pid << " started new worker # " << nthreads << endl;
1833 }
1834 steal_ctr++;
1835 for(int j = 0; j < 10; j++) {
1836 sleep(1);
1837 }
1838 }
1839 }
1840 #endif
1841
1842 #if (__cplusplus >= 201103L)
1843 while(all_threads_done < nthreads) {
1844 SLEEP(10);
1845 }
1846 #else
1847 for (int i = 0; i < nthreads - 1; i++) {
1848 threads[i]->join();
1849 }
1850 #endif
1851
1852 #ifndef _WIN32
1853 if(thread_stealing) {
1854 del_pid(thread_stealing_dir.c_str(), pid);
1855 }
1856 #endif
1857 }
1858 if(refs != NULL) delete refs;
1859
1860 for (int i = 0; i < nthreads - 1; i++) {
1861 if (threads[i] != NULL) {
1862 delete threads[i];
1863 }
1864 }
1865 }
1866
1867 #define SWITCH_TO_FW_INDEX() { \
1868 /* Evict the mirror index from memory if necessary */ \
1869 if(ebwtBw.isInMemory()) ebwtBw.evictFromMemory(); \
1870 assert(!ebwtBw.isInMemory()); \
1871 /* Load the forward index into memory if necessary */ \
1872 if(!ebwtFw.isInMemory()) { \
1873 Timer _t(cerr, "Time loading forward index: ", timing); \
1874 ebwtFw.loadIntoMemory(-1, !noRefNames, startVerbose); \
1875 } \
1876 assert(ebwtFw.isInMemory()); \
1877 _patsrc.reset(); /* rewind pattern source to first pattern */ \
1878 }
1879
1880 #define SWITCH_TO_BW_INDEX() { \
1881 /* Evict the forward index from memory if necessary */ \
1882 if(ebwtFw.isInMemory()) ebwtFw.evictFromMemory(); \
1883 assert(!ebwtFw.isInMemory()); \
1884 /* Load the forward index into memory if necessary */ \
1885 if(!ebwtBw.isInMemory()) { \
1886 Timer _t(cerr, "Time loading mirror index: ", timing); \
1887 ebwtBw.loadIntoMemory(-1, !noRefNames, startVerbose); \
1888 } \
1889 assert(ebwtBw.isInMemory()); \
1890 _patsrc.reset(); /* rewind pattern source to first pattern */ \
1891 }
1892
1893 #define ASSERT_NO_HITS_FW(ebwtfw) \
1894 if(sanityCheck && os.size() > 0) { \
1895 EList<Hit> hits; \
1896 uint32_t threeRevOff = (seedMms <= 3) ? s : 0; \
1897 uint32_t twoRevOff = (seedMms <= 2) ? s : 0; \
1898 uint32_t oneRevOff = (seedMms <= 1) ? s : 0; \
1899 uint32_t unrevOff = (seedMms == 0) ? s : 0; \
1900 if(hits.size() > 0) { \
1901 /* Print offending hit obtained by oracle */ \
1902 ::printHit( \
1903 os, \
1904 hits[0], \
1905 patFw, \
1906 plen, \
1907 unrevOff, \
1908 oneRevOff, \
1909 twoRevOff, \
1910 threeRevOff, \
1911 ebwtfw); /* ebwtFw */ \
1912 } \
1913 assert_eq(0, hits.size()); \
1914 }
1915
1916 #define ASSERT_NO_HITS_RC(ebwtfw) \
1917 if(sanityCheck && os.size() > 0) { \
1918 EList<Hit> hits; \
1919 uint32_t threeRevOff = (seedMms <= 3) ? s : 0; \
1920 uint32_t twoRevOff = (seedMms <= 2) ? s : 0; \
1921 uint32_t oneRevOff = (seedMms <= 1) ? s : 0; \
1922 uint32_t unrevOff = (seedMms == 0) ? s : 0; \
1923 if(hits.size() > 0) { \
1924 /* Print offending hit obtained by oracle */ \
1925 ::printHit( \
1926 os, \
1927 hits[0], \
1928 patRc, \
1929 plen, \
1930 unrevOff, \
1931 oneRevOff, \
1932 twoRevOff, \
1933 threeRevOff, \
1934 ebwtfw); /* ebwtFw */ \
1935 } \
1936 assert_eq(0, hits.size()); \
1937 }
1938
1939 static PatternComposer* twoOrThreeMismatchSearch_patsrc;
1940 static HitSink* twoOrThreeMismatchSearch_sink;
1941 static Ebwt* twoOrThreeMismatchSearch_ebwtFw;
1942 static Ebwt* twoOrThreeMismatchSearch_ebwtBw;
1943 static EList<BTRefString >* twoOrThreeMismatchSearch_os;
1944 static SyncBitset* twoOrThreeMismatchSearch_doneMask;
1945 static SyncBitset* twoOrThreeMismatchSearch_hitMask;
1946 static bool twoOrThreeMismatchSearch_two;
1947 static BitPairReference* twoOrThreeMismatchSearch_refs;
1948
1949
1950 /**
1951 * A statefulness-aware worker driver. Uses UnpairedExactAlignerV1.
1952 */
1953 #if (__cplusplus >= 201103L)
1954 //void twoOrThreeMismatchSearchWorkerStateful::operator()() const {
1955 static void twoOrThreeMismatchSearchWorkerStateful(void *vp) {
1956 thread_tracking_pair *p = (thread_tracking_pair*) vp;
1957 int tid = p->tid;
1958 #else
1959 static void twoOrThreeMismatchSearchWorkerStateful(void *vp) {
1960 int tid = *((int*)vp);
1961 #endif
1962 if(thread_stealing) {
1963 increment_thread_counter();
1964 }
1965 PatternComposer& _patsrc = *twoOrThreeMismatchSearch_patsrc;
1966 HitSink& _sink = *twoOrThreeMismatchSearch_sink;
1967 Ebwt& ebwtFw = *twoOrThreeMismatchSearch_ebwtFw;
1968 Ebwt& ebwtBw = *twoOrThreeMismatchSearch_ebwtBw;
1969 EList<BTRefString >& os = *twoOrThreeMismatchSearch_os;
1970 BitPairReference* refs = twoOrThreeMismatchSearch_refs;
1971 static bool two = twoOrThreeMismatchSearch_two;
1972
1973 // Global initialization
1974 PatternSourcePerThreadFactory* patsrcFact = createPatsrcFactory(_patsrc, tid, readsPerBatch);
1975 HitSinkPerThreadFactory* sinkFact = createSinkFactory(_sink, tid);
1976
1977 ChunkPool *pool = new ChunkPool(chunkSz * 1024, chunkPoolMegabytes * 1024 * 1024, chunkVerbose);
1978 Unpaired23mmAlignerV1Factory alSEfact(
1979 ebwtFw,
1980 &ebwtBw,
1981 two,
1982 !nofw,
1983 !norc,
1984 _sink,
1985 *sinkFact,
1986 NULL, //&cacheFw,
1987 NULL, //&cacheBw,
1988 cacheLimit,
1989 pool,
1990 refs,
1991 os,
1992 !noMaqRound,
1993 !better,
1994 strandFix,
1995 rangeMode,
1996 verbose,
1997 quiet,
1998 seed);
1999 Paired23mmAlignerV1Factory alPEfact(
2000 ebwtFw,
2001 &ebwtBw,
2002 !nofw,
2003 !norc,
2004 useV1,
2005 two,
2006 _sink,
2007 *sinkFact,
2008 mate1fw,
2009 mate2fw,
2010 minInsert,
2011 maxInsert,
2012 dontReconcileMates,
2013 mhits, // for symCeiling
2014 mixedThresh,
2015 mixedAttemptLim,
2016 NULL, //&cacheFw,
2017 NULL, //&cacheBw,
2018 cacheLimit,
2019 pool,
2020 refs, os,
2021 reportSe,
2022 !noMaqRound,
2023 !better,
2024 strandFix,
2025 rangeMode,
2026 verbose,
2027 quiet,
2028 seed);
2029 {
2030 MixedMultiAligner multi(
2031 prefetchWidth,
2032 qUpto,
2033 alSEfact,
2034 alPEfact,
2035 *patsrcFact);
2036 // Run that mother
2037 multi.run(false, tid);
2038 // MultiAligner must be destroyed before patsrcFact
2039 }
2040
2041 #if (__cplusplus >= 201103L)
2042 p->done->fetch_add(1);
2043 #endif
2044 if(thread_stealing) {
2045 decrement_thread_counter();
2046 }
2047
2048 delete patsrcFact;
2049 delete sinkFact;
2050 delete pool;
2051 return;
2052 }
2053
2054 #if (__cplusplus >= 201103L)
2055 //void twoOrThreeMismatchSearchWorkerFull::operator()() const {
2056 static void twoOrThreeMismatchSearchWorkerFull(void *vp) {
2057 thread_tracking_pair *p = (thread_tracking_pair*) vp;
2058 int tid = p->tid;
2059 #else
2060 static void twoOrThreeMismatchSearchWorkerFull(void *vp) {
2061 int tid = *((int*)vp);
2062 #endif
2063 if(thread_stealing) {
2064 increment_thread_counter();
2065 }
2066 PatternComposer& _patsrc = *twoOrThreeMismatchSearch_patsrc;
2067 HitSink& _sink = *twoOrThreeMismatchSearch_sink;
2068 EList<BTRefString >& os = *twoOrThreeMismatchSearch_os;
2069 bool two = twoOrThreeMismatchSearch_two;
2070 PatternSourcePerThreadFactory* patsrcFact = createPatsrcFactory(_patsrc, tid, readsPerBatch);
2071 PatternSourcePerThread* patsrc = patsrcFact->create();
2072 HitSinkPerThreadFactory* sinkFact = createSinkFactory(_sink, tid);
2073 HitSinkPerThread* sink = sinkFact->create();
2074 /* Per-thread initialization */
2075 EbwtSearchParams params(
2076 *sink, /* HitSink */
2077 os, /* reference sequences */
2078 true, /* read is forward */
2079 true); /* index is forward */
2080 Ebwt& ebwtFw = *twoOrThreeMismatchSearch_ebwtFw;
2081 Ebwt& ebwtBw = *twoOrThreeMismatchSearch_ebwtBw;
2082 GreedyDFSRangeSource btr1(
2083 &ebwtFw, params,
2084 0xffffffff, // qualThresh
2085 // Do not impose maximums in 2/3-mismatch mode
2086 0xffffffff, // max backtracks (no limit)
2087 0, // reportPartials (don't)
2088 true, // reportExacts
2089 rangeMode, // reportRanges
2090 NULL, // seedlings
2091 NULL, // mutations
2092 verbose, // verbose
2093 &os,
2094 false); // considerQuals
2095 GreedyDFSRangeSource bt2(
2096 &ebwtBw, params,
2097 0xffffffff, // qualThresh
2098 // Do not impose maximums in 2/3-mismatch mode
2099 0xffffffff, // max backtracks (no limit)
2100 0, // reportPartials (no)
2101 true, // reportExacts
2102 rangeMode, // reportRanges
2103 NULL, // seedlings
2104 NULL, // mutations
2105 verbose, // verbose
2106 &os,
2107 false); // considerQuals
2108 GreedyDFSRangeSource bt3(
2109 &ebwtFw, params,
2110 0xffffffff, // qualThresh (none)
2111 // Do not impose maximums in 2/3-mismatch mode
2112 0xffffffff, // max backtracks (no limit)
2113 0, // reportPartials (don't)
2114 true, // reportExacts
2115 rangeMode, // reportRanges
2116 NULL, // seedlings
2117 NULL, // mutations
2118 verbose, // verbose
2119 &os,
2120 false); // considerQuals
2121 GreedyDFSRangeSource bthh3(
2122 &ebwtFw, params,
2123 0xffffffff, // qualThresh
2124 // Do not impose maximums in 2/3-mismatch mode
2125 0xffffffff, // max backtracks (no limit)
2126 0, // reportPartials (don't)
2127 true, // reportExacts
2128 rangeMode, // reportRanges
2129 NULL, // seedlings
2130 NULL, // mutations
2131 verbose, // verbose
2132 &os,
2133 false, // considerQuals
2134 true); // halfAndHalf
2135 bool skipped = false;
2136 pair<bool, bool> get_read_ret = make_pair(false, false);
2137 #ifdef PER_THREAD_TIMING
2138 uint64_t ncpu_changeovers = 0;
2139 uint64_t nnuma_changeovers = 0;
2140
2141 int current_cpu = 0, current_node = 0;
2142 get_cpu_and_node(current_cpu, current_node);
2143
2144 std::stringstream ss;
2145 std::string msg;
2146 ss << "thread: " << tid << " time: ";
2147 msg = ss.str();
2148 {
2149 Timer timer(std::cout, msg.c_str());
2150 #endif
2151 while(true) { // Read read-in loop
2152 #ifdef PER_THREAD_TIMING
2153 int cpu = 0, node = 0;
2154 get_cpu_and_node(cpu, node);
2155 if(cpu != current_cpu) {
2156 ncpu_changeovers++;
2157 current_cpu = cpu;
2158 }
2159 if(node != current_node) {
2160 nnuma_changeovers++;
2161 current_node = node;
2162 }
2163 #endif
2164 FINISH_READ(patsrc);
2165 GET_READ(patsrc);
2166 patid += 0; // kill unused variable warning
2167 uint32_t plen = (uint32_t)patFw.length();
2168 uint32_t s = plen;
2169 uint32_t s3 = s >> 1; // length of 3' half of seed
2170 uint32_t s5 = (s >> 1) + (s & 1); // length of 5' half of seed
2171 #define DONEMASK_SET(p)
2172 #include "search_23mm_phase1.c"
2173 #include "search_23mm_phase2.c"
2174 #include "search_23mm_phase3.c"
2175 #undef DONEMASK_SET
2176 }
2177 FINISH_READ(patsrc);
2178 if(thread_stealing) {
2179 decrement_thread_counter();
2180 }
2181 #ifdef PER_THREAD_TIMING
2182 ss.str("");
2183 ss.clear();
2184 ss << "thread: " << tid << " cpu_changeovers: " << ncpu_changeovers << std::endl
2185 << "thread: " << tid << " node_changeovers: " << nnuma_changeovers << std::endl;
2186 std::cout << ss.str();
2187 }
2188 #endif
2189 #if (__cplusplus >= 201103L)
2190 p->done->fetch_add(1);
2191 #endif
2192 // Threads join at end of Phase 1
2193 WORKER_EXIT();
2194 }
2195
2196 static void twoOrThreeMismatchSearchFull(
2197 PatternComposer& _patsrc, /// pattern source
2198 HitSink& _sink, /// hit sink
2199 Ebwt& ebwtFw, /// index of original text
2200 Ebwt& ebwtBw, /// index of mirror text
2201 EList<BTRefString >& os, /// text strings, if available (empty otherwise)
2202 bool two = true) /// true -> 2, false -> 3
2203 {
2204 // Global initialization
2205 assert(!ebwtFw.isInMemory());
2206 assert(!ebwtBw.isInMemory());
2207 {
2208 // Load the other half of the index into memory
2209 Timer _t(cerr, "Time loading forward index: ", timing);
2210 ebwtFw.loadIntoMemory(-1, !noRefNames, startVerbose);
2211 }
2212 {
2213 // Load the other half of the index into memory
2214 Timer _t(cerr, "Time loading mirror index: ", timing);
2215 ebwtBw.loadIntoMemory(-1, !noRefNames, startVerbose);
2216 }
2217 // Create range caches, which are shared among all aligners
2218 BitPairReference *refs = NULL;
2219 bool pair = mates1.size() > 0 || mates12.size() > 0;
2220 if(pair && mixedThresh < 0xffffffff) {
2221 Timer _t(cerr, "Time loading reference: ", timing);
2222 refs = new BitPairReference(adjustedEbwtFileBase, sanityCheck, NULL, &os, false, true, useMm, useShmem, mmSweep, verbose, startVerbose);
2223 if(!refs->loaded()) throw 1;
2224 }
2225 twoOrThreeMismatchSearch_refs = refs;
2226 twoOrThreeMismatchSearch_patsrc = &_patsrc;
2227 twoOrThreeMismatchSearch_sink = &_sink;
2228 twoOrThreeMismatchSearch_ebwtFw = &ebwtFw;
2229 twoOrThreeMismatchSearch_ebwtBw = &ebwtBw;
2230 twoOrThreeMismatchSearch_os = &os;
2231 twoOrThreeMismatchSearch_doneMask = NULL;
2232 twoOrThreeMismatchSearch_hitMask = NULL;
2233 twoOrThreeMismatchSearch_two = two;
2234
2235 int tids[max(nthreads, thread_ceiling)];
2236 #if (__cplusplus >= 201103L)
2237 EList<std::thread*> threads;
2238 EList<thread_tracking_pair> tps;
2239 tps.resizeExact(max(nthreads, thread_ceiling));
2240 #else
2241 EList<tthread::thread*> threads;
2242 #endif
2243
2244 #if (__cplusplus >= 201103L)
2245 std::atomic<int> all_threads_done;
2246 all_threads_done = 0;
2247 #endif
2248
2249 CHUD_START();
2250 {
2251 Timer _t(cerr, "End-to-end 2/3-mismatch full-index search: ", timing);
2252
2253 #ifndef _WIN32
2254 int pid = 0;
2255 if(thread_stealing) {
2256 pid = getpid();
2257 write_pid(thread_stealing_dir.c_str(), pid);
2258 thread_counter = 0;
2259 }
2260 #endif
2261
2262 for(int i = 0; i < nthreads; i++) {
2263 tids[i] = i;
2264 #if (__cplusplus >= 201103L)
2265 tps[i].tid = tids[i];
2266 tps[i].done = &all_threads_done;
2267 if (i == nthreads - 1) {
2268 if(stateful) {
2269 twoOrThreeMismatchSearchWorkerStateful((void*)&tps[i]);
2270 } else {
2271 twoOrThreeMismatchSearchWorkerFull((void*)&tps[i]);
2272 }
2273 }
2274 else {
2275 if(stateful) {
2276 threads.push_back(new std::thread(twoOrThreeMismatchSearchWorkerStateful, (void*)&tps[i]));
2277 } else {
2278 threads.push_back(new std::thread(twoOrThreeMismatchSearchWorkerFull, (void*)&tps[i]));
2279 }
2280 threads[i]->detach();
2281 SLEEP(10);
2282 }
2283 #else
2284 if (i == nthreads - 1) {
2285 if(stateful) {
2286 twoOrThreeMismatchSearchWorkerStateful((void*)(tids + i));
2287 } else {
2288 twoOrThreeMismatchSearchWorkerFull((void*)(tids + i));
2289 }
2290 }
2291 else {
2292 if(stateful) {
2293 threads.push_back(new tthread::thread(twoOrThreeMismatchSearchWorkerStateful, (void *)(tids + i)));
2294 } else {
2295 threads.push_back(new tthread::thread(twoOrThreeMismatchSearchWorkerFull, (void *)(tids + i)));
2296 }
2297 }
2298 #endif
2299 }
2300
2301 #ifndef _WIN32
2302 if(thread_stealing) {
2303 int orig_threads = nthreads, steal_ctr = 1;
2304 for(int j = 0; j < 10; j++) {
2305 sleep(1);
2306 }
2307 while(thread_counter > 0) {
2308 if(steal_thread(pid, orig_threads)) {
2309 nthreads++;
2310 tids[nthreads-1] = nthreads;
2311 #if (__cplusplus >= 201103L)
2312 tps[nthreads-1].tid = nthreads - 1;
2313 tps[nthreads-1].done = &all_threads_done;
2314 if(stateful) {
2315 threads.push_back(new std::thread(twoOrThreeMismatchSearchWorkerStateful, (void*)&tps[nthreads-1]));
2316 } else {
2317 threads.push_back(new std::thread(twoOrThreeMismatchSearchWorkerFull, (void*)&tps[nthreads-1]));
2318 }
2319 threads[nthreads-1]->detach();
2320 SLEEP(10);
2321 #else
2322 if(stateful) {
2323 threads.push_back(new tthread::thread(twoOrThreeMismatchSearchWorkerStateful, (void *)(tids + nthreads - 1)));
2324 } else {
2325 threads.push_back(new tthread::thread(twoOrThreeMismatchSearchWorkerFull, (void *)(tids + nthreads - 1)));
2326 }
2327 #endif
2328 cerr << "pid " << pid << " started new worker # " << nthreads << endl;
2329 }
2330 steal_ctr++;
2331 for(int j = 0; j < 10; j++) {
2332 sleep(1);
2333 }
2334 }
2335 }
2336 #endif
2337
2338 #if (__cplusplus >= 201103L)
2339 while(all_threads_done < nthreads) {
2340 SLEEP(10);
2341 }
2342 #else
2343 for (int i = 0; i < nthreads - 1; i++) {
2344 threads[i]->join();
2345 }
2346 #endif
2347
2348 #ifndef _WIN32
2349 if(thread_stealing) {
2350 del_pid(thread_stealing_dir.c_str(), pid);
2351 }
2352 #endif
2353 }
2354 if(refs != NULL) delete refs;
2355
2356 for (int i = 0; i < nthreads - 1; i++) {
2357 if (threads[i] != NULL) {
2358 delete threads[i];
2359 }
2360 }
2361 return;
2362 }
2363
2364 static PatternComposer* seededQualSearch_patsrc;
2365 static HitSink* seededQualSearch_sink;
2366 static Ebwt* seededQualSearch_ebwtFw;
2367 static Ebwt* seededQualSearch_ebwtBw;
2368 static EList<BTRefString >* seededQualSearch_os;
2369 static SyncBitset* seededQualSearch_doneMask;
2370 static SyncBitset* seededQualSearch_hitMask;
2371 static PartialAlignmentManager* seededQualSearch_pamFw;
2372 static PartialAlignmentManager* seededQualSearch_pamRc;
2373 static int seededQualSearch_qualCutoff;
2374 static BitPairReference* seededQualSearch_refs;
2375
2376 #if (__cplusplus >= 201103L)
2377 //void seededQualSearchWorkerFull::operator()() const {
2378 static void seededQualSearchWorkerFull(void *vp) {
2379 thread_tracking_pair *p = (thread_tracking_pair*) vp;
2380 int tid = p->tid;
2381 #else
2382 static void seededQualSearchWorkerFull(void *vp) {
2383 int tid = *((int*)vp);
2384 #endif
2385 if(thread_stealing) {
2386 increment_thread_counter();
2387 }
2388 PatternComposer& _patsrc = *seededQualSearch_patsrc;
2389 HitSink& _sink = *seededQualSearch_sink;
2390 EList<BTRefString >& os = *seededQualSearch_os;
2391 int qualCutoff = seededQualSearch_qualCutoff;
2392 PatternSourcePerThreadFactory* patsrcFact = createPatsrcFactory(_patsrc, tid, readsPerBatch);
2393 PatternSourcePerThread* patsrc = patsrcFact->create();
2394 HitSinkPerThreadFactory* sinkFact = createSinkFactory(_sink, tid);
2395 HitSinkPerThread* sink = sinkFact->create();
2396 /* Per-thread initialization */
2397 EbwtSearchParams params(
2398 *sink, /* HitSink */
2399 os, /* reference sequences */
2400 true, /* read is forward */
2401 true); /* index is forward */
2402 Ebwt& ebwtFw = *seededQualSearch_ebwtFw;
2403 Ebwt& ebwtBw = *seededQualSearch_ebwtBw;
2404 PartialAlignmentManager * pamRc = NULL;
2405 PartialAlignmentManager * pamFw = NULL;
2406 if(seedMms > 0) {
2407 pamRc = new PartialAlignmentManager(64);
2408 pamFw = new PartialAlignmentManager(64);
2409 }
2410 EList<PartialAlignment> pals;
2411 // GreedyDFSRangeSource for finding exact hits for the forward-
2412 // oriented read
2413 GreedyDFSRangeSource btf1(
2414 &ebwtFw, params,
2415 qualCutoff, // qualThresh
2416 maxBtsBetter, // max backtracks
2417 0, // reportPartials (don't)
2418 true, // reportExacts
2419 rangeMode, // reportRanges
2420 NULL, // seedlings
2421 NULL, // mutations
2422 verbose, // verbose
2423 &os,
2424 false); // considerQuals
2425 GreedyDFSRangeSource bt1(
2426 &ebwtFw, params,
2427 qualCutoff, // qualThresh
2428 maxBtsBetter, // max backtracks
2429 0, // reportPartials (don't)
2430 true, // reportExacts
2431 rangeMode, // reportRanges
2432 NULL, // seedlings
2433 NULL, // mutations
2434 verbose, // verbose
2435 &os, // reference sequences
2436 true, // considerQuals
2437 false, !noMaqRound);
2438 // GreedyDFSRangeSource to search for hits for cases 1F, 2F, 3F
2439 GreedyDFSRangeSource btf2(
2440 &ebwtBw, params,
2441 qualCutoff, // qualThresh
2442 maxBtsBetter, // max backtracks
2443 0, // reportPartials (no)
2444 true, // reportExacts
2445 rangeMode, // reportRanges
2446 NULL, // partial alignment manager
2447 NULL, // mutations
2448 verbose, // verbose
2449 &os, // reference sequences
2450 true, // considerQuals
2451 false, !noMaqRound);
2452 // GreedyDFSRangeSource to search for partial alignments for case 4R
2453 GreedyDFSRangeSource btr2(
2454 &ebwtBw, params,
2455 qualCutoff, // qualThresh (none)
2456 maxBtsBetter, // max backtracks
2457 seedMms, // report partials (up to seedMms mms)
2458 true, // reportExacts
2459 rangeMode, // reportRanges
2460 pamRc, // partial alignment manager
2461 NULL, // mutations
2462 verbose, // verbose
2463 &os, // reference sequences
2464 true, // considerQuals
2465 false, !noMaqRound);
2466 // GreedyDFSRangeSource to search for seedlings for case 4F
2467 GreedyDFSRangeSource btf3(
2468 &ebwtFw, params,
2469 qualCutoff, // qualThresh (none)
2470 maxBtsBetter, // max backtracks
2471 seedMms, // reportPartials (do)
2472 true, // reportExacts
2473 rangeMode, // reportRanges
2474 pamFw, // seedlings
2475 NULL, // mutations
2476 verbose, // verbose
2477 &os, // reference sequences
2478 true, // considerQuals
2479 false, !noMaqRound);
2480 // GreedyDFSRangeSource to search for hits for case 4R by extending
2481 // the partial alignments found in Phase 2
2482 GreedyDFSRangeSource btr3(
2483 &ebwtFw, params,
2484 qualCutoff, // qualThresh
2485 maxBtsBetter, // max backtracks
2486 0, // reportPartials (don't)
2487 true, // reportExacts
2488 rangeMode,// reportRanges
2489 NULL, // seedlings
2490 NULL, // mutations
2491 verbose, // verbose
2492 &os, // reference sequences
2493 true, // considerQuals
2494 false, !noMaqRound);
2495 // The half-and-half GreedyDFSRangeSource
2496 GreedyDFSRangeSource btr23(
2497 &ebwtFw, params,
2498 qualCutoff, // qualThresh
2499 maxBtsBetter, // max backtracks
2500 0, // reportPartials (don't)
2501 true, // reportExacts
2502 rangeMode,// reportRanges
2503 NULL, // seedlings
2504 NULL, // mutations
2505 verbose, // verbose
2506 &os,
2507 true, // considerQuals
2508 true, // halfAndHalf
2509 !noMaqRound);
2510 // GreedyDFSRangeSource to search for hits for case 4F by extending
2511 // the partial alignments found in Phase 3
2512 GreedyDFSRangeSource btf4(
2513 &ebwtBw, params,
2514 qualCutoff, // qualThresh
2515 maxBtsBetter, // max backtracks
2516 0, // reportPartials (don't)
2517 true, // reportExacts
2518 rangeMode,// reportRanges
2519 NULL, // seedlings
2520 NULL, // mutations
2521 verbose, // verbose
2522 &os, // reference sequences
2523 true, // considerQuals
2524 false, !noMaqRound);
2525 // Half-and-half GreedyDFSRangeSource for forward read
2526 GreedyDFSRangeSource btf24(
2527 &ebwtBw, params,
2528 qualCutoff, // qualThresh
2529 maxBtsBetter, // max backtracks
2530 0, // reportPartials (don't)
2531 true, // reportExacts
2532 rangeMode,// reportRanges
2533 NULL, // seedlings
2534 NULL, // mutations
2535 verbose, // verbose
2536 &os,
2537 true, // considerQuals
2538 true, // halfAndHalf
2539 !noMaqRound);
2540 EList<QueryMutation> muts;
2541 bool skipped = false;
2542 pair<bool, bool> get_read_ret = make_pair(false, false);
2543 #ifdef PER_THREAD_TIMING
2544 uint64_t ncpu_changeovers = 0;
2545 uint64_t nnuma_changeovers = 0;
2546
2547 int current_cpu = 0, current_node = 0;
2548 get_cpu_and_node(current_cpu, current_node);
2549
2550 std::stringstream ss;
2551 std::string msg;
2552 ss << "thread: " << tid << " time: ";
2553 msg = ss.str();
2554 {
2555 Timer timer(std::cout, msg.c_str());
2556 #endif
2557 while(true) {
2558 #ifdef PER_THREAD_TIMING
2559 int cpu = 0, node = 0;
2560 get_cpu_and_node(cpu, node);
2561 if(cpu != current_cpu) {
2562 ncpu_changeovers++;
2563 current_cpu = cpu;
2564 }
2565 if(node != current_node) {
2566 nnuma_changeovers++;
2567 current_node = node;
2568 }
2569 #endif
2570 FINISH_READ(patsrc);
2571 GET_READ(patsrc);
2572 uint32_t plen = (uint32_t)patFw.length();
2573 uint32_t s = seedLen;
2574 uint32_t s3 = (s >> 1); /* length of 3' half of seed */
2575 uint32_t s5 = (s >> 1) + (s & 1); /* length of 5' half of seed */
2576 uint32_t qs = min<uint32_t>(plen, s);
2577 uint32_t qs3 = qs >> 1;
2578 uint32_t qs5 = (qs >> 1) + (qs & 1);
2579 #define DONEMASK_SET(p)
2580 #include "search_seeded_phase1.c"
2581 #include "search_seeded_phase2.c"
2582 #include "search_seeded_phase3.c"
2583 #include "search_seeded_phase4.c"
2584 #undef DONEMASK_SET
2585 }
2586 FINISH_READ(patsrc);
2587 if(thread_stealing) {
2588 decrement_thread_counter();
2589 }
2590 if(seedMms > 0) {
2591 delete pamRc;
2592 delete pamFw;
2593 }
2594 #ifdef PER_THREAD_TIMING
2595 ss.str("");
2596 ss.clear();
2597 ss << "thread: " << tid << " cpu_changeovers: " << ncpu_changeovers << std::endl
2598 << "thread: " << tid << " node_changeovers: " << nnuma_changeovers << std::endl;
2599 std::cout << ss.str();
2600 }
2601 #endif
2602 #if (__cplusplus >= 201103L)
2603 p->done->fetch_add(1);
2604 #endif
2605 WORKER_EXIT();
2606 }
2607 #if (__cplusplus >= 201103L)
2608 //void seededQualSearchWorkerFullStateful::operator()() const {
2609 static void seededQualSearchWorkerFullStateful(void *vp) {
2610 thread_tracking_pair *p = (thread_tracking_pair*) vp;
2611 int tid = p->tid;
2612 #else
2613 static void seededQualSearchWorkerFullStateful(void *vp) {
2614 int tid = *((int*)vp);
2615 #endif
2616 if(thread_stealing) {
2617 increment_thread_counter();
2618 }
2619 PatternComposer& _patsrc = *seededQualSearch_patsrc;
2620 HitSink& _sink = *seededQualSearch_sink;
2621 Ebwt& ebwtFw = *seededQualSearch_ebwtFw;
2622 Ebwt& ebwtBw = *seededQualSearch_ebwtBw;
2623 EList<BTRefString >& os = *seededQualSearch_os;
2624 int qualCutoff = seededQualSearch_qualCutoff;
2625 BitPairReference* refs = seededQualSearch_refs;
2626
2627 // Global initialization
2628 PatternSourcePerThreadFactory* patsrcFact = createPatsrcFactory(_patsrc, tid, readsPerBatch);
2629 HitSinkPerThreadFactory* sinkFact = createSinkFactory(_sink, tid);
2630 ChunkPool *pool = new ChunkPool(chunkSz * 1024, chunkPoolMegabytes * 1024 * 1024, chunkVerbose);
2631
2632 AlignerMetrics *metrics = NULL;
2633 if(stats) {
2634 metrics = new AlignerMetrics();
2635 }
2636 UnpairedSeedAlignerFactory alSEfact(
2637 ebwtFw,
2638 &ebwtBw,
2639 !nofw,
2640 !norc,
2641 seedMms,
2642 seedLen,
2643 qualCutoff,
2644 maxBts,
2645 _sink,
2646 *sinkFact,
2647 NULL, //&cacheFw,
2648 NULL, //&cacheBw,
2649 cacheLimit,
2650 pool,
2651 refs,
2652 os,
2653 !noMaqRound,
2654 !better,
2655 strandFix,
2656 rangeMode,
2657 verbose,
2658 quiet,
2659 seed,
2660 metrics);
2661 PairedSeedAlignerFactory alPEfact(
2662 ebwtFw,
2663 &ebwtBw,
2664 useV1,
2665 !nofw,
2666 !norc,
2667 seedMms,
2668 seedLen,
2669 qualCutoff,
2670 maxBts,
2671 _sink,
2672 *sinkFact,
2673 mate1fw,
2674 mate2fw,
2675 minInsert,
2676 maxInsert,
2677 dontReconcileMates,
2678 mhits, // for symCeiling
2679 mixedThresh,
2680 mixedAttemptLim,
2681 NULL, //&cacheFw,
2682 NULL, //&cacheBw,
2683 cacheLimit,
2684 pool,
2685 refs,
2686 os,
2687 reportSe,
2688 !noMaqRound,
2689 !better,
2690 strandFix,
2691 rangeMode,
2692 verbose,
2693 quiet,
2694 seed);
2695 {
2696 MixedMultiAligner multi(
2697 prefetchWidth,
2698 qUpto,
2699 alSEfact,
2700 alPEfact,
2701 *patsrcFact);
2702 // Run that mother
2703 multi.run(false, tid);
2704 // MultiAligner must be destroyed before patsrcFact
2705 }
2706 if(metrics != NULL) {
2707 metrics->printSummary();
2708 delete metrics;
2709 }
2710
2711 #if (__cplusplus >= 201103L)
2712 p->done->fetch_add(1);
2713 #endif
2714 if(thread_stealing) {
2715 decrement_thread_counter();
2716 }
2717
2718 delete patsrcFact;
2719 delete sinkFact;
2720 delete pool;
2721 return;
2722 }
2723
2724 /**
2725 * Search for a good alignments for each read using criteria that
2726 * correspond somewhat faithfully to Maq's. Search is aided by a pair
2727 * of Ebwt indexes, one for the original references, and one for the
2728 * transpose of the references. Neither index should be loaded upon
2729 * entry to this function.
2730 *
2731 * Like Maq, we treat the first 24 base pairs of the read (those
2732 * closest to the 5' end) differently from the remainder of the read.
2733 * We call the first 24 base pairs the "seed."
2734 */
2735 static void seededQualCutoffSearchFull(
2736 int seedLen, /// length of seed (not a maq option)
2737 int qualCutoff, /// maximum sum of mismatch qualities
2738 /// like maq map's -e option
2739 /// default: 70
2740 int seedMms, /// max # mismatches allowed in seed
2741 /// (like maq map's -n option)
2742 /// Can only be 1 or 2, default: 1
2743 PatternComposer& _patsrc, /// pattern source
2744 HitSink& _sink, /// hit sink
2745 Ebwt& ebwtFw, /// index of original text
2746 Ebwt& ebwtBw, /// index of mirror text
2747 EList<BTRefString >& os) /// text strings, if available (empty otherwise)
2748 {
2749 // Global intialization
2750 assert_leq(seedMms, 3);
2751
2752 seededQualSearch_patsrc = &_patsrc;
2753 seededQualSearch_sink = &_sink;
2754 seededQualSearch_ebwtFw = &ebwtFw;
2755 seededQualSearch_ebwtBw = &ebwtBw;
2756 seededQualSearch_os = &os;
2757 seededQualSearch_doneMask = NULL;
2758 seededQualSearch_hitMask = NULL;
2759 seededQualSearch_pamFw = NULL;
2760 seededQualSearch_pamRc = NULL;
2761 seededQualSearch_qualCutoff = qualCutoff;
2762
2763 // Create range caches, which are shared among all aligners
2764 BitPairReference *refs = NULL;
2765 bool pair = mates1.size() > 0 || mates12.size() > 0;
2766 if(pair && mixedThresh < 0xffffffff) {
2767 Timer _t(cerr, "Time loading reference: ", timing);
2768 refs = new BitPairReference(adjustedEbwtFileBase, sanityCheck, NULL, &os, false, true, useMm, useShmem, mmSweep, verbose, startVerbose);
2769 if(!refs->loaded()) throw 1;
2770 }
2771 seededQualSearch_refs = refs;
2772
2773 int tids[max(nthreads, thread_ceiling)];
2774 #if (__cplusplus >= 201103L)
2775 EList<std::thread*> threads;
2776 EList<thread_tracking_pair> tps;
2777 tps.resizeExact(max(nthreads, thread_ceiling));
2778 #else
2779 EList<tthread::thread*> threads;
2780 #endif
2781
2782 #if (__cplusplus >= 201103L)
2783 std::atomic<int> all_threads_done;
2784 all_threads_done = 0;
2785 #endif
2786
2787 SWITCH_TO_FW_INDEX();
2788 assert(!ebwtBw.isInMemory());
2789 {
2790 // Load the other half of the index into memory
2791 Timer _t(cerr, "Time loading mirror index: ", timing);
2792 ebwtBw.loadIntoMemory(-1, !noRefNames, startVerbose);
2793 }
2794 CHUD_START();
2795 {
2796 // Phase 1: Consider cases 1R and 2R
2797 Timer _t(cerr, "Seeded quality full-index search: ", timing);
2798
2799 #ifndef _WIN32
2800 int pid = 0;
2801 if(thread_stealing) {
2802 pid = getpid();
2803 write_pid(thread_stealing_dir.c_str(), pid);
2804 thread_counter = 0;
2805 }
2806 #endif
2807
2808 for(int i = 0; i < nthreads; i++) {
2809 tids[i] = i;
2810 #if (__cplusplus >= 201103L)
2811 tps[i].tid = tids[i];
2812 tps[i].done = &all_threads_done;
2813 if (i == nthreads - 1) {
2814 if(stateful) {
2815 seededQualSearchWorkerFullStateful((void*)&tps[i]);
2816 } else {
2817 seededQualSearchWorkerFull((void*)&tps[i]);
2818 }
2819 }
2820 else {
2821 if(stateful) {
2822 threads.push_back(new std::thread(seededQualSearchWorkerFullStateful, (void*)&tps[i]));
2823 } else {
2824 threads.push_back(new std::thread(seededQualSearchWorkerFull, (void*)&tps[i]));
2825 }
2826 threads[i]->detach();
2827 SLEEP(10);
2828 }
2829 #else
2830 if (i == nthreads - 1) {
2831 if(stateful) {
2832 seededQualSearchWorkerFullStateful((void*)(tids + i));
2833 } else {
2834 seededQualSearchWorkerFull((void*)(tids + i));
2835 }
2836 }
2837 else {
2838 if(stateful) {
2839 threads.push_back(new tthread::thread(seededQualSearchWorkerFullStateful, (void *)(tids + i)));
2840 } else {
2841 threads.push_back(new tthread::thread(seededQualSearchWorkerFull, (void *)(tids + i)));
2842 }
2843 }
2844 #endif
2845 }
2846
2847 #ifndef _WIN32
2848 if(thread_stealing) {
2849 int orig_threads = nthreads, steal_ctr = 1;
2850 for(int j = 0; j < 10; j++) {
2851 sleep(1);
2852 }
2853 while(thread_counter > 0) {
2854 if(steal_thread(pid, orig_threads)) {
2855 nthreads++;
2856 tids[nthreads-1] = nthreads - 1;
2857 #if (__cplusplus >= 201103L)
2858 tps[nthreads-1].tid = nthreads - 1;
2859 tps[nthreads-1].done = &all_threads_done;
2860 if(stateful) {
2861 threads.push_back(new std::thread(seededQualSearchWorkerFullStateful, (void*)&tps[nthreads-1]));
2862 } else {
2863 threads.push_back(new std::thread(seededQualSearchWorkerFull, (void*)&tps[nthreads-1]));
2864 }
2865 threads[nthreads-1]->detach();
2866 SLEEP(10);
2867 #else
2868 if(stateful) {
2869 threads.push_back(new tthread::thread(seededQualSearchWorkerFullStateful, (void *)(tids + nthreads - 1)));
2870 } else {
2871 threads.push_back(new tthread::thread(seededQualSearchWorkerFull, (void *)(tids + nthreads - 1)));
2872 }
2873 #endif
2874 cerr << "pid " << pid << " started new worker # " << nthreads << endl;
2875 }
2876 steal_ctr++;
2877 for(int j = 0; j < 10; j++) {
2878 sleep(1);
2879 }
2880 }
2881 }
2882 #endif
2883
2884 #if (__cplusplus >= 201103L)
2885 while(all_threads_done < nthreads) {
2886 SLEEP(10);
2887 }
2888 #else
2889 for (int i = 0; i < nthreads - 1; i++) {
2890 threads[i]->join();
2891 }
2892 #endif
2893
2894 #ifndef _WIN32
2895 if(thread_stealing) {
2896 del_pid(thread_stealing_dir.c_str(), pid);
2897 }
2898 #endif
2899 }
2900
2901 if(refs != NULL) {
2902 delete refs;
2903 }
2904
2905 for (int i = 0; i < nthreads - 1; i++) {
2906 if (threads[i] != NULL) {
2907 delete threads[i];
2908 }
2909 }
2910
2911 ebwtBw.evictFromMemory();
2912 }
2913
2914 /**
2915 * Return a new dynamically allocated PatternSource for the given
2916 * format, using the given list of strings as the filenames to read
2917 * from or as the sequences themselves (i.e. if -c was used).
2918 */
2919 static PatternSource*
2920 patsrcFromStrings(int format,
2921 const EList<string>& reads,
2922 const EList<string>* quals)
2923 {
2924 switch(format) {
2925 case FASTA:
2926 return new FastaPatternSource (reads, quals,
2927 trim3, trim5);
2928 case FASTA_CONT:
2929 return new FastaContinuousPatternSource (
2930 reads, fastaContLen,
2931 fastaContFreq);
2932 case RAW:
2933 return new RawPatternSource (reads, trim3, trim5);
2934 case FASTQ:
2935 return new FastqPatternSource (reads, trim3, trim5,
2936 solexaQuals, phred64Quals,
2937 integerQuals);
2938 case INTERLEAVED:
2939 return new FastqPatternSource (reads, trim3, trim5,
2940 solexaQuals, phred64Quals,
2941 integerQuals, true /* is interleaved */);
2942 case TAB_MATE:
2943 return new TabbedPatternSource(reads, false, trim3, trim5);
2944 case CMDLINE:
2945 return new VectorPatternSource(reads, trim3, trim5);
2946 default: {
2947 cerr << "Internal error; bad patsrc format: " << format << endl;
2948 throw 1;
2949 }
2950 }
2951 }
2952
2953 static string argstr;
2954
2955 static void driver(const char * type,
2956 const string& ebwtFileBase,
2957 const string& query,
2958 const EList<string>& queries,
2959 const EList<string>& qualities,
2960 const string& outfile)
2961 {
2962 if(verbose || startVerbose) {
2963 cerr << "Entered driver(): "; logTime(cerr, true);
2964 }
2965 // Vector of the reference sequences; used for sanity-checking
2966 EList<BTRefString> os;
2967 // Read reference sequences from the command-line or from a FASTA file
2968 if(!origString.empty()) {
2969 // Determine if it's a file by looking at whether it has a FASTA-like
2970 // extension
2971 size_t len = origString.length();
2972 if((len >= 6 && origString.substr(len-6) == ".fasta") ||
2973 (len >= 4 && origString.substr(len-4) == ".mfa") ||
2974 (len >= 4 && origString.substr(len-4) == ".fas") ||
2975 (len >= 4 && origString.substr(len-4) == ".fna") ||
2976 (len >= 3 && origString.substr(len-3) == ".fa"))
2977 {
2978 // Read fasta file
2979 EList<string> origFiles;
2980 tokenize(origString, ",", origFiles);
2981 readSequenceFiles(origFiles, os);
2982 } else {
2983 // Read sequence
2984 readSequenceString(origString, os);
2985 }
2986 }
2987 // Adjust
2988 adjustedEbwtFileBase = adjustEbwtBase(argv0, ebwtFileBase, verbose);
2989 bool isBt2Index = false;
2990 if (gEbwt_ext == "bt2" || gEbwt_ext == "bt2l") {
2991 isBt2Index = true;
2992 }
2993
2994
2995 EList<PatternSource*> patsrcs_a;
2996 EList<PatternSource*> patsrcs_b;
2997 EList<PatternSource*> patsrcs_ab;
2998
2999 // If there were any first-mates specified, we will operate in
3000 // stateful mode
3001 bool paired = mates1.size() > 0 || mates12.size() > 0;
3002 if(paired) stateful = true;
3003
3004 // Create list of pattern sources for paired reads appearing
3005 // interleaved in a single file
3006 if(verbose || startVerbose) {
3007 cerr << "Creating paired-end patsrcs: "; logTime(cerr, true);
3008 }
3009 for(size_t i = 0; i < mates12.size(); i++) {
3010 const EList<string>* qs = &mates12;
3011 EList<string> tmp;
3012 if(fileParallel) {
3013 // Feed query files one to each PatternSource
3014 qs = &tmp;
3015 tmp.push_back(mates12[i]);
3016 assert_eq(1, tmp.size());
3017 }
3018 patsrcs_ab.push_back(patsrcFromStrings(format, *qs, NULL));
3019 if(!fileParallel) {
3020 break;
3021 }
3022 }
3023
3024 // Create list of pattern sources for paired reads
3025 for(size_t i = 0; i < mates1.size(); i++) {
3026 const EList<string>* qs = &mates1;
3027 const EList<string>* quals = &qualities1;
3028 EList<string> tmpSeq;
3029 EList<string> tmpQual;
3030 if(fileParallel) {
3031 // Feed query files one to each PatternSource
3032 qs = &tmpSeq;
3033 tmpSeq.push_back(mates1[i]);
3034 quals = &tmpSeq;
3035 tmpQual.push_back(qualities1[i]);
3036 assert_eq(1, tmpSeq.size());
3037 }
3038 if(quals->empty()) quals = NULL;
3039 patsrcs_a.push_back(patsrcFromStrings(format, *qs, quals));
3040 if(!fileParallel) {
3041 break;
3042 }
3043 }
3044
3045 // Create list of pattern sources for paired reads
3046 for(size_t i = 0; i < mates2.size(); i++) {
3047 const EList<string>* qs = &mates2;
3048 const EList<string>* quals = &qualities2;
3049 EList<string> tmpSeq;
3050 EList<string> tmpQual;
3051 if(fileParallel) {
3052 // Feed query files one to each PatternSource
3053 qs = &tmpSeq;
3054 tmpSeq.push_back(mates2[i]);
3055 quals = &tmpQual;
3056 tmpQual.push_back(qualities2[i]);
3057 assert_eq(1, tmpSeq.size());
3058 }
3059 if(quals->empty()) quals = NULL;
3060 patsrcs_b.push_back(patsrcFromStrings(format, *qs, quals));
3061 if(!fileParallel) {
3062 break;
3063 }
3064 }
3065 // All mates/mate files must be paired
3066 assert_eq(patsrcs_a.size(), patsrcs_b.size());
3067
3068 // Create list of pattern sources for the unpaired reads
3069 if(verbose || startVerbose) {
3070 cerr << "Creating single-end patsrcs: "; logTime(cerr, true);
3071 }
3072 for(size_t i = 0; i < queries.size(); i++) {
3073 const EList<string>* qs = &queries;
3074 const EList<string>* quals = &qualities;
3075 PatternSource* patsrc = NULL;
3076 EList<string> tmpSeq;
3077 EList<string> tmpQual;
3078 if(fileParallel) {
3079 // Feed query files one to each PatternSource
3080 qs = &tmpSeq;
3081 tmpSeq.push_back(queries[i]);
3082 quals = &tmpQual;
3083 tmpQual.push_back(qualities[i]);
3084 assert_eq(1, tmpSeq.size());
3085 }
3086 if(quals->empty()) quals = NULL;
3087 patsrc = patsrcFromStrings(format, *qs, quals);
3088 assert(patsrc != NULL);
3089 patsrcs_a.push_back(patsrc);
3090 patsrcs_b.push_back(NULL);
3091 if(!fileParallel) {
3092 break;
3093 }
3094 }
3095
3096 if(verbose || startVerbose) {
3097 cerr << "Creating PatternSource: "; logTime(cerr, true);
3098 }
3099 PatternComposer *patsrc = NULL;
3100 if(mates12.size() > 0) {
3101 patsrc = new SoloPatternComposer(patsrcs_ab);
3102 } else {
3103 patsrc = new DualPatternComposer(patsrcs_a, patsrcs_b);
3104 }
3105
3106 // Open hit output file
3107 if(verbose || startVerbose) {
3108 cerr << "Opening hit output file: "; logTime(cerr, true);
3109 }
3110 OutFileBuf *fout;
3111 if(!outfile.empty()) {
3112 fout = new OutFileBuf(outfile.c_str(), false);
3113 } else {
3114 fout = new OutFileBuf();
3115 }
3116 // Initialize Ebwt object and read in header
3117 if(verbose || startVerbose) {
3118 cerr << "About to initialize fw Ebwt: "; logTime(cerr, true);
3119 }
3120 Ebwt ebwt(adjustedEbwtFileBase,
3121 -1, // don't care about entireReverse
3122 true, // index is for the forward direction
3123 /* overriding: */ offRate,
3124 /* overriding: */ isaRate,
3125 useMm, // whether to use memory-mapped files
3126 useShmem, // whether to use shared memory
3127 mmSweep, // sweep memory-mapped files
3128 !noRefNames, // load names?
3129 verbose, // whether to be talkative
3130 startVerbose, // talkative during initialization
3131 false /*passMemExc*/,
3132 sanityCheck,
3133 isBt2Index);
3134 Ebwt* ebwtBw = NULL;
3135 // We need the mirror index if mismatches are allowed
3136 if(mismatches > 0 || maqLike) {
3137 if(verbose || startVerbose) {
3138 cerr << "About to initialize rev Ebwt: "; logTime(cerr, true);
3139 }
3140 ebwtBw = new Ebwt(
3141 adjustedEbwtFileBase + ".rev",
3142 -1, // don't care about entireReverse
3143 false, // index is for the reverse direction
3144 /* overriding: */ offRate,
3145 /* overriding: */ isaRate,
3146 useMm, // whether to use memory-mapped files
3147 useShmem, // whether to use shared memory
3148 mmSweep, // sweep memory-mapped files
3149 !noRefNames, // load names?
3150 verbose, // whether to be talkative
3151 startVerbose, // talkative during initialization
3152 false /*passMemExc*/,
3153 sanityCheck,
3154 isBt2Index);
3155 }
3156 if(!os.empty()) {
3157 for(size_t i = 0; i < os.size(); i++) {
3158 size_t olen = os[i].length();
3159 int longestStretch = 0;
3160 int curStretch = 0;
3161 for(size_t j = 0; j < olen; j++) {
3162 if((int)os[i][j] < 4) {
3163 curStretch++;
3164 if(curStretch > longestStretch) longestStretch = curStretch;
3165 } else {
3166 curStretch = 0;
3167 }
3168 }
3169 if(longestStretch < 1) {
3170 os.erase(i--);
3171 }
3172 }
3173 }
3174 if(sanityCheck && !os.empty()) {
3175 // Sanity check number of patterns and pattern lengths in Ebwt
3176 // against original strings
3177 assert_eq(os.size(), ebwt.nPat());
3178 for(size_t i = 0; i < os.size(); i++) {
3179 assert_eq(os[i].length(), ebwt.plen()[i]);
3180 }
3181 ebwt.loadIntoMemory(-1, !noRefNames, startVerbose);
3182 ebwt.checkOrigs(os, false);
3183 ebwt.evictFromMemory();
3184 }
3185 {
3186 Timer _t(cerr, "Time searching: ", timing);
3187 if(verbose || startVerbose) {
3188 cerr << "Creating HitSink: "; logTime(cerr, true);
3189 }
3190 // Set up hit sink; if sanityCheck && !os.empty() is true,
3191 // then instruct the sink to "retain" hits in a vector in
3192 // memory so that we can easily sanity check them later on
3193 HitSink *sink;
3194 EList<string>* refnames = &ebwt.refnames();
3195 if(noRefNames) refnames = NULL;
3196 if(outType == OUTPUT_FULL) {
3197 sink = new VerboseHitSink(
3198 *fout, offBase,
3199 printCost,
3200 suppressOuts,
3201 fullRef,
3202 dumpAlBase,
3203 dumpUnalBase,
3204 dumpMaxBase,
3205 format == TAB_MATE, sampleMax,
3206 refnames, nthreads,
3207 outBatchSz, partitionSz);
3208 } else if(outType == OUTPUT_SAM) {
3209 SAMHitSink *sam = new SAMHitSink(
3210 *fout,
3211 fullRef, samNoQnameTrunc,
3212 dumpAlBase,
3213 dumpUnalBase,
3214 dumpMaxBase,
3215 format == TAB_MATE,
3216 sampleMax,
3217 refnames,
3218 nthreads,
3219 outBatchSz,
3220 reorder);
3221 if(!samNoHead) {
3222 EList<string> refnames;
3223 if(!samNoSQ) {
3224 readEbwtRefnames(adjustedEbwtFileBase, refnames);
3225 }
3226 sam->appendHeaders(
3227 sam->out(),
3228 ebwt.nPat(),
3229 refnames, samNoSQ,
3230 ebwt.plen(), fullRef,
3231 samNoQnameTrunc,
3232 argstr.c_str(),
3233 rgs.empty() ? NULL : rgs.c_str());
3234 }
3235 sink = sam;
3236 } else {
3237 cerr << "Invalid output type: " << outType << endl;
3238 throw 1;
3239 }
3240 if(verbose || startVerbose) {
3241 cerr << "Dispatching to search driver: "; logTime(cerr, true);
3242 }
3243 if(maqLike) {
3244 seededQualCutoffSearchFull(seedLen,
3245 qualThresh,
3246 seedMms,
3247 *patsrc,
3248 *sink,
3249 ebwt, // forward index
3250 *ebwtBw, // mirror index (not optional)
3251 os); // references, if available
3252 }
3253 else if(mismatches > 0) {
3254 if(mismatches == 1) {
3255 assert(ebwtBw != NULL);
3256 mismatchSearchFull(*patsrc, *sink, ebwt, *ebwtBw, os);
3257 } else if(mismatches == 2 || mismatches == 3) {
3258 twoOrThreeMismatchSearchFull(*patsrc, *sink, ebwt, *ebwtBw, os, mismatches == 2);
3259 } else {
3260 cerr << "Error: " << mismatches << " is not a supported number of mismatches" << endl;
3261 throw 1;
3262 }
3263 } else {
3264 // Search without mismatches
3265 // Note that --fast doesn't make a difference here because
3266 // we're only loading half of the index anyway
3267 exactSearch(*patsrc, *sink, ebwt, os);
3268 }
3269 // Evict any loaded indexes from memory
3270 if(ebwt.isInMemory()) {
3271 ebwt.evictFromMemory();
3272 }
3273 if(ebwtBw != NULL) {
3274 delete ebwtBw;
3275 }
3276 sink->finish(hadoopOut); // end the hits section of the hit file
3277 for(size_t i = 0; i < patsrcs_a.size(); i++) {
3278 assert(patsrcs_a[i] != NULL);
3279 delete patsrcs_a[i];
3280 }
3281 for(size_t i = 0; i < patsrcs_b.size(); i++) {
3282 if(patsrcs_b[i] != NULL) {
3283 delete patsrcs_b[i];
3284 }
3285 }
3286 for(size_t i = 0; i < patsrcs_ab.size(); i++) {
3287 if(patsrcs_ab[i] != NULL) {
3288 delete patsrcs_ab[i];
3289 }
3290 }
3291 delete patsrc;
3292 delete sink;
3293 if(fout != NULL) delete fout;
3294 }
3295 }
3296
3297 // C++ name mangling is disabled for the bowtie() function to make it
3298 // easier to use Bowtie as a library.
3299 extern "C" {
3300
3301 /**
3302 * Main bowtie entry function. Parses argc/argv style command-line
3303 * options, sets global configuration variables, and calls the driver()
3304 * function.
3305 */
3306 int bowtie(int argc, const char **argv) {
3307 try {
3308 #if (__cplusplus >= 201103L)
3309 #ifdef WITH_AFFINITY
3310 //CWILKS: adjust this depending on # of hyperthreads per core
3311 pinning_observer pinner( 2 /* the number of hyper threads on each core */ );
3312 pinner.observe( true );
3313 #endif
3314 #endif
3315 // Reset all global state, including getopt state
3316 opterr = optind = 1;
3317 resetOptions();
3318 for(int i = 0; i < argc; i++) {
3319 argstr += argv[i];
3320 if(i < argc-1) argstr += " ";
3321 }
3322 string query; // read query string(s) from this file
3323 EList<string> queries;
3324 string outfile; // write query results to this file
3325 if(startVerbose) { cerr << "Entered main(): "; logTime(cerr, true); }
3326 parseOptions(argc, argv);
3327 argv0 = argv[0];
3328 if(showVersion) {
3329 cout << argv0 << " version " << BOWTIE_VERSION << endl;
3330 if(sizeof(void*) == 4) {
3331 cout << "32-bit" << endl;
3332 } else if(sizeof(void*) == 8) {
3333 cout << "64-bit" << endl;
3334 } else {
3335 cout << "Neither 32- nor 64-bit: sizeof(void*) = " << sizeof(void*) << endl;
3336 }
3337 cout << "Built on " << BUILD_HOST << endl;
3338 cout << BUILD_TIME << endl;
3339 cout << "Compiler: " << COMPILER_VERSION << endl;
3340 cout << "Options: " << COMPILER_OPTIONS << endl;
3341 cout << "Sizeof {int, long, long long, void*, size_t, off_t}: {"
3342 << sizeof(int)
3343 << ", " << sizeof(long) << ", " << sizeof(long long)
3344 << ", " << sizeof(void *) << ", " << sizeof(size_t)
3345 << ", " << sizeof(off_t) << "}" << endl;
3346 return 0;
3347 }
3348 #ifdef CHUD_PROFILING
3349 chudInitialize();
3350 chudAcquireRemoteAccess();
3351 #endif
3352 {
3353 Timer _t(cerr, "Overall time: ", timing);
3354 if(startVerbose) {
3355 cerr << "Parsing index and read arguments: "; logTime(cerr, true);
3356 }
3357
3358 // Get index basename
3359 if (ebwtFile.empty()) {
3360 if(optind >= argc) {
3361 cerr << "No index, query, or output file specified!" << endl;
3362 printUsage(cerr);
3363 return 1;
3364 }
3365 std::cerr << "Setting the index via positional argument"
3366 << " will be deprecated in a future release."
3367 << " Please use -x option instead." << std::endl;
3368 ebwtFile = argv[optind++];
3369 }
3370
3371 // Get query filename
3372 if(optind >= argc) {
3373 if(mates1.size() > 0 || mates12.size() > 0) {
3374 query = "";
3375 } else {
3376 cerr << "No query or output file specified!" << endl;
3377 printUsage(cerr);
3378 return 1;
3379 }
3380 } else if (mates1.size() == 0 && mates12.size() == 0) {
3381 query = argv[optind++];
3382 // Tokenize the list of query files
3383 tokenize(query, ",", queries);
3384 if(queries.size() < 1) {
3385 cerr << "Tokenized query file list was empty!" << endl;
3386 printUsage(cerr);
3387 return 1;
3388 }
3389 }
3390
3391 // Get output filename
3392 if(optind < argc) {
3393 outfile = argv[optind++];
3394 }
3395
3396 // Extra parametesr?
3397 if(optind < argc) {
3398 cerr << "Extra parameter(s) specified: ";
3399 for(int i = optind; i < argc; i++) {
3400 cerr << "\"" << argv[i] << "\"";
3401 if(i < argc-1) cerr << ", ";
3402 }
3403 cerr << endl;
3404 if(mates1.size() > 0) {
3405 cerr << "Note that if <mates> files are specified using -1/-2, a <singles> file cannot" << endl
3406 << "also be specified. Please run bowtie separately for mates and singles." << endl;
3407 }
3408 throw 1;
3409 }
3410
3411 // Optionally summarize
3412 if(verbose) {
3413 cout << "Input ebwt file: \"" << ebwtFile << "\"" << endl;
3414 cout << "Query inputs (DNA, " << file_format_names[format] << "):" << endl;
3415 for(size_t i = 0; i < queries.size(); i++) {
3416 cout << " " << queries[i] << endl;
3417 }
3418 cout << "Quality inputs:" << endl;
3419 for(size_t i = 0; i < qualities.size(); i++) {
3420 cout << " " << qualities[i] << endl;
3421 }
3422 cout << "Output file: \"" << outfile << "\"" << endl;
3423 cout << "Local endianness: " << (currentlyBigEndian()? "big":"little") << endl;
3424 cout << "Sanity checking: " << (sanityCheck? "enabled":"disabled") << endl;
3425 #ifdef NDEBUG
3426 cout << "Assertions: disabled" << endl;
3427 #else
3428 cout << "Assertions: enabled" << endl;
3429 #endif
3430 }
3431 if(ipause) {
3432 cout << "Press key to continue..." << endl;
3433 getchar();
3434 }
3435 driver("DNA", ebwtFile, query, queries, qualities, outfile);
3436 CHUD_STOP();
3437 }
3438 #ifdef CHUD_PROFILING
3439 chudReleaseRemoteAccess();
3440 #endif
3441 #if (__cplusplus >= 201103L)
3442 #ifdef WITH_AFFINITY
3443 // Always disable observation before observers destruction
3444 //tracker.observe( false );
3445 pinner.observe( false );
3446 #endif
3447 #endif
3448 return 0;
3449 } catch(exception& e) {
3450 cerr << "Command: ";
3451 for(int i = 0; i < argc; i++) cerr << argv[i] << " ";
3452 cerr << endl;
3453 return 1;
3454 } catch(int e) {
3455 if(e != 0) {
3456 cerr << "Command: ";
3457 for(int i = 0; i < argc; i++) cerr << argv[i] << " ";
3458 cerr << endl;
3459 }
3460 return e;
3461 }
3462 } // bowtie()
3463 } // extern "C"
3464