Imported Upstream version 0.12.7
[bowtie.git] / ref_aligner.h
1 /*
2  * ref_aligner.h
3  */
4
5 #ifndef REF_ALIGNER_H_
6 #define REF_ALIGNER_H_
7
8 #include <stdint.h>
9 #include <iostream>
10 #include <vector>
11 #include <stdexcept>
12 #include "seqan/sequence.h"
13 #include "alphabet.h"
14 #include "range.h"
15 #include "reference.h"
16
17 // Let the reference-aligner buffer size be 16K by default.  If more
18 // room is required, a new buffer must be allocated from the heap.
19 const static int REF_ALIGNER_BUFSZ = 16 * 1024;
20
21 /**
22  * Abstract parent class for classes that look for alignments by
23  * matching against the reference sequence directly.  This is useful
24  * both for sanity-checking results from the Bowtie index and for
25  * finding mates when the reference location of the opposite mate is
26  * known.
27  */
28 template<typename TStr>
29 class RefAligner {
30
31         typedef seqan::String<seqan::Dna5> TDna5Str;
32         typedef seqan::String<char> TCharStr;
33         typedef std::vector<uint32_t> TU32Vec;
34         typedef std::vector<Range> TRangeVec;
35         typedef std::pair<uint64_t, uint64_t> TU64Pair;
36         typedef std::set<TU64Pair> TSetPairs;
37
38 public:
39         RefAligner(bool color,
40                    bool verbose = false,
41                    bool quiet = false,
42                    uint32_t seedLen = 0,
43                    uint32_t qualMax = 0xffffffff,
44                    bool maqPenalty = false) :
45                 color_(color), verbose_(verbose), seedLen_(seedLen),
46                 qualMax_(qualMax), maqPenalty_(maqPenalty), refbuf_(buf_),
47                 refbufSz_(REF_ALIGNER_BUFSZ), freeRefbuf_(false)
48                 { }
49
50         /**
51          * Free the reference-space alignment buffer if this object
52          * allocated it.
53          */
54         virtual ~RefAligner() {
55                 if(freeRefbuf_) {
56                         delete[] refbuf_;
57                 }
58         }
59
60         /**
61          * Find one alignment of qry:quals in the range begin-end in
62          * reference string ref.  Store the alignment details in range.
63          */
64         virtual void find(uint32_t numToFind,
65                           const uint32_t tidx,
66                           const BitPairReference *refs,
67                           const TDna5Str& qry,
68                           const TCharStr& quals,
69                           uint32_t begin,
70                           uint32_t end,
71                           TRangeVec& ranges,
72                           TU32Vec& results,
73                           TSetPairs* pairs = NULL,
74                           uint32_t aoff = 0xffffffff,
75                           bool seedOnLeft = false)
76         {
77                 assert_gt(numToFind, 0);
78                 assert_gt(end, begin);
79                 uint32_t spread = end - begin + (color_ ? 1 : 0);
80                 uint32_t spreadPlus = spread + 12;
81                 // Make sure the buffer is large enough to accommodate the spread
82                 if(spreadPlus > this->refbufSz_) {
83                         this->newBuf(spreadPlus);
84                 }
85                 // Read in the relevant stretch of the reference string
86                 int offset = refs->getStretch(this->refbuf_, tidx, begin, spread);
87                 uint8_t *buf = ((uint8_t*)this->refbuf_) + offset;
88                 if(color_) {
89                         // Colorize buffer
90                         for(size_t i = 0; i < (end-begin); i++) {
91                                 assert_leq((int)buf[i], 4);
92                                 buf[i] = dinuc2color[(int)buf[i]][(int)buf[i+1]];
93                         }
94                 }
95                 // Look for alignments
96                 ASSERT_ONLY(uint32_t irsz = ranges.size());
97                 anchor64Find(numToFind, tidx, buf, qry, quals, begin,
98                              end, ranges, results, pairs, aoff, seedOnLeft);
99 #ifndef NDEBUG
100                 for(size_t i = irsz; i < results.size(); i++) {
101                         assert_eq(ranges[i].numMms, ranges[i].mms.size());
102                         assert_eq(ranges[i].numMms, ranges[i].refcs.size());
103                 }
104 #endif
105         }
106
107         /**
108          * Find one alignment of qry:quals in the range begin-end in
109          * reference string ref.  Store the alignment details in range.
110          * Uses a combination of the anchor bases and 64-bit arithmetic to
111          * find anchors quickly.
112          */
113         virtual void anchor64Find(uint32_t numToFind,
114                         uint32_t tidx,
115                         uint8_t* ref,
116                     const TDna5Str& qry,
117                     const TCharStr& quals,
118                     uint32_t begin,
119                     uint32_t end,
120                         TRangeVec& ranges,
121                         TU32Vec& results,
122                         TSetPairs* pairs = NULL,
123                         uint32_t aoff = 0xffffffff,
124                         bool seedOnLeft = false) const = 0;
125
126         /**
127          * Set a new reference-sequence buffer.
128          */
129         void setBuf(uint32_t *newbuf, uint32_t newsz) {
130                 if(freeRefbuf_) {
131                         delete[] refbuf_;
132                         freeRefbuf_ = false;
133                 }
134                 refbuf_ = newbuf;
135                 refbufSz_ = newsz;
136         }
137
138         /**
139          * Set a new reference-sequence buffer.
140          */
141         void newBuf(uint32_t newsz) {
142                 if(freeRefbuf_) {
143                         delete[] refbuf_;
144                 }
145                 try {
146                         refbuf_ = new uint32_t[(newsz + 3) / 4];
147                         if(refbuf_ == NULL) throw std::bad_alloc();
148                 } catch(std::bad_alloc& e) {
149                         cerr << "Error: Could not allocate reference-space alignment buffer of " << newsz << "B" << endl;
150                         throw 1;
151                 }
152                 refbufSz_ = newsz;
153                 freeRefbuf_ = true;
154         }
155
156 protected:
157         bool      color_;     /// whether to colorize reference buffers before handing off
158         bool      verbose_;   /// be talkative
159         uint32_t  seedLen_;   /// length of seed region for read
160         uint32_t  qualMax_;   /// maximum sum of quality penalties
161         bool      maqPenalty_;/// whether to round as in Maq
162         uint32_t *refbuf_;    /// pointer to current reference buffer
163         uint32_t  refbufSz_;  /// size of current reference buffer
164         uint32_t  buf_[REF_ALIGNER_BUFSZ / 4]; /// built-in reference buffer (may be superseded)
165         bool      freeRefbuf_; /// whether refbuf_ points to something we should delete
166 };
167
168 /**
169  * Concrete RefAligner for finding nearby exact hits given an anchor
170  * hit.
171  */
172 template<typename TStr>
173 class ExactRefAligner : public RefAligner<TStr> {
174
175         typedef seqan::String<seqan::Dna5> TDna5Str;
176         typedef seqan::String<char> TCharStr;
177         typedef std::vector<uint32_t> TU32Vec;
178         typedef std::vector<Range> TRangeVec;
179         typedef std::pair<uint64_t, uint64_t> TU64Pair;
180         typedef std::set<TU64Pair> TSetPairs;
181
182 public:
183
184         ExactRefAligner(bool color, bool verbose, bool quiet) :
185                 RefAligner<TStr>(color, verbose, quiet) { }
186
187         virtual ~ExactRefAligner() { }
188
189 protected:
190         /**
191          * Because we're doing end-to-end exact, we don't care which end of
192          * 'qry' is the 5' end.
193          */
194         void naiveFind(uint32_t numToFind,
195                        uint32_t tidx,
196                        uint8_t* ref,
197                    const TDna5Str& qry,
198                    const TCharStr& quals,
199                    uint32_t begin,
200                    uint32_t end,
201                        TRangeVec& ranges,
202                        TU32Vec& results,
203                        TSetPairs* pairs,
204                        uint32_t aoff,
205                        bool seedOnLeft) const
206         {
207                 assert_gt(numToFind, 0);
208                 const uint32_t qlen = seqan::length(qry);
209                 assert_geq(end - begin, qlen); // caller should have checked this
210                 assert_gt(end, begin);
211                 assert_gt(qlen, 0);
212                 uint32_t qend = end - qlen;
213                 uint32_t lim = qend - begin;
214                 uint32_t halfway = begin + (lim >> 1);
215                 bool hi = false;
216                 for(uint32_t i = 1; i <= lim+1; i++) {
217                         uint32_t ri;  // leftmost position in candidate alignment
218                         uint32_t rir; // same, minus begin; for indexing into ref[]
219                         if(hi) {
220                                 ri = halfway + (i >> 1); rir = ri - begin;
221                                 assert_leq(ri, qend);
222                         } else {
223                                 ri = halfway - (i >> 1); rir = ri - begin;
224                                 assert_geq(ri, begin);
225                         }
226                         hi = !hi;
227                         // Do the naive comparison
228                         bool match = true;
229                         for(uint32_t j = 0; j < qlen; j++) {
230 #if 0
231                                 // Count Ns in the reference as mismatches
232                                 const int q = (int)qry[j];
233                                 const int r = (int)ref[rir + j];
234                                 assert_leq(q, 4);
235                                 assert_leq(r, 4);
236                                 if(q == 4 || r == 4 || q != r) {
237                                         // Mismatch
238                                         match = false;
239                                         break;
240                                 }
241                                 // Match; continue
242 #else
243                                 // Disallow alignments that involve an N in the
244                                 // reference
245                                 const int r = (int)ref[rir + j];
246                                 if(r & 4) {
247                                         // N in reference; bail
248                                         match = false;
249                                         break;
250                                 }
251                                 const int q = (int)qry[j];
252                                 assert_leq(q, 4);
253                                 assert_lt(r, 4);
254                                 if(q != r) {
255                                         // Mismatch
256                                         match = false;
257                                         break;
258                                 }
259                                 // Match; continue
260 #endif
261                         }
262                         if(match) {
263                                 ranges.resize(ranges.size()+1);
264                                 Range& range = ranges.back();
265                                 range.stratum = 0;
266                                 range.numMms = 0;
267                                 assert_eq(0, range.mms.size());
268                                 assert_eq(0, range.refcs.size());
269                                 results.push_back(ri);
270                         }
271                 }
272                 return; // no match
273         }
274
275         /**
276          * Because we're doing end-to-end exact, we don't care which end of
277          * 'qry' is the 5' end.
278          */
279         virtual void anchor64Find(uint32_t numToFind,
280                     uint32_t tidx,
281                         uint8_t *ref,
282                     const TDna5Str& qry,
283                     const TCharStr& quals,
284                     uint32_t begin,
285                     uint32_t end,
286                         TRangeVec& ranges,
287                         TU32Vec& results,
288                         TSetPairs* pairs,
289                         uint32_t aoff, // offset of anchor mate
290                         bool seedOnLeft) const
291         {
292                 assert_gt(numToFind, 0);
293                 ASSERT_ONLY(const uint32_t rangesInitSz = ranges.size());
294                 ASSERT_ONLY(uint32_t duplicates = 0);
295                 ASSERT_ONLY(uint32_t r2i = 0);
296                 const uint32_t qlen = seqan::length(qry);
297                 assert_geq(end - begin, qlen); // caller should have checked this
298                 assert_gt(end, begin);
299                 assert_gt(qlen, 0);
300 #ifndef NDEBUG
301                 // Get all naive hits
302                 TRangeVec r2; TU32Vec re2;
303                 naiveFind(numToFind, tidx, ref, qry, quals, begin, end, r2,
304                           re2, pairs, aoff, seedOnLeft);
305 #endif
306                 const uint32_t anchorBitPairs = min<int>(qlen, 32);
307                 // anchorOverhang = # read bases not included in the anchor
308                 const uint32_t anchorOverhang = qlen <= 32 ? 0 : qlen - 32;
309                 const uint32_t lim = end - qlen - begin;
310                 const uint32_t halfway = begin + (lim >> 1);
311                 uint64_t anchor = 0llu;
312                 uint64_t buffw = 0llu;
313                 // Set up a mask that we'll apply to the two bufs every round
314                 // to discard bits that were rotated out of the anchor area
315                 uint64_t clearMask = 0xffffffffffffffffllu;
316                 bool useMask = false;
317                 if(anchorBitPairs < 32) {
318                         clearMask >>= ((32-anchorBitPairs) << 1);
319                         useMask = true;
320                 }
321                 const int lhsShift = ((anchorBitPairs - 1) << 1);
322                 // Build the contents of the 'anchor' dword and the initial
323                 // contents of the 'buffw' dword.  If there are fewer than 32
324                 // anchorBitPairs, the content will be packed into the least
325                 // significant bits of the word.
326                 uint32_t skipLeftToRights = 0;
327                 uint32_t skipRightToLefts = 0;
328                 for(uint32_t i = 0; i < anchorBitPairs; i++) {
329                         int c = (int)qry[i]; // next query character
330                         assert_leq(c, 4);
331                         if(c & 4) {
332                                 assert_eq(r2.size(), ranges.size() - rangesInitSz);
333                                 return; // can't match if query has Ns
334                         }
335                         int r = (int)ref[halfway - begin + i]; // next reference character
336                         if(r & 4) {
337                                 r = 0;
338                                 // The right-to-left direction absorbs the candidate
339                                 // alignment based at 'halfway'; so skipLeftToRights is
340                                 // i, not i+1
341                                 skipLeftToRights = max(skipLeftToRights, i);
342                                 skipRightToLefts = max(skipRightToLefts, anchorBitPairs - i);
343                         }
344                         assert_lt(r, 4);
345                         assert_lt(c, 4);
346                         anchor  = ((anchor  << 2llu) | c);
347                         buffw = ((buffw << 2llu) | r);
348                 }
349                 // Check whether read is disqualified by Ns outside of the anchor
350                 // region
351                 for(uint32_t i = anchorBitPairs; i < qlen; i++) {
352                         if((int)qry[i] == 4) {
353                                 assert_eq(r2.size(), ranges.size() - rangesInitSz);
354                                 return; // can't match if query has Ns
355                         }
356                 }
357                 uint64_t bufbw = buffw;
358                 // We're moving the right-hand edge of the anchor along until
359                 // it's 'anchorOverhang' chars from the end of the target region.
360                 // Note that we're not making a 3'/5' distinction here; if we
361                 // were, we might need to make the 'anchorOverhang' adjustment on
362                 // the left end of the range rather than the right.
363                 bool hi = false;
364                 uint32_t riHi  = halfway;
365                 uint32_t rirHi = halfway - begin;
366                 uint32_t rirHiAnchor = rirHi + anchorBitPairs - 1;
367                 uint32_t riLo  = halfway + 1;
368                 uint32_t rirLo = halfway - begin + 1;
369                 for(uint32_t i = 1; i <= lim + 1; i++) {
370                         int r;       // new reference char
371                         assert_lt(skipLeftToRights, qlen);
372                         assert_leq(skipRightToLefts, qlen);
373                         if(hi) {
374                                 hi = false;
375                                 // Moving left-to-right
376                                 riHi++; rirHi++; rirHiAnchor++;
377                                 r = (int)ref[rirHiAnchor];
378                                 if(r & 4) {
379                                         r = 0;
380                                         skipLeftToRights = anchorBitPairs;
381                                 }
382                                 assert_lt(r, 4);
383                                 // Bring in new base pair at the least significant
384                                 // position
385                                 buffw = ((buffw << 2llu) | r);
386                                 if(useMask) buffw &= clearMask;
387                                 if(skipLeftToRights > 0) {
388                                         skipLeftToRights--;
389                                         continue;
390                                 }
391                                 if(buffw != anchor) {
392                                         continue;
393                                 }
394                         } else {
395                                 hi = true;
396                                 // Moving right-to-left
397                                 riLo--; rirLo--;
398                                 r = (int)ref[rirLo];
399                                 if(r & 4) {
400                                         r = 0;
401                                         skipRightToLefts = qlen;
402                                 }
403                                 assert_lt(r, 4);
404                                 if(i >= 2) {
405                                         bufbw >>= 2llu;
406                                         // Bring in new base pair at the most significant
407                                         // position
408                                         bufbw |= ((uint64_t)r << lhsShift);
409                                 }
410                                 if(skipRightToLefts > 0) {
411                                         skipRightToLefts--;
412                                         continue;
413                                 }
414                                 if(bufbw != anchor) {
415                                         continue;
416                                 }
417                         }
418                         // Seed hit!
419                         bool foundHit = true;
420                         uint32_t ri = hi ? riLo : riHi;
421                         uint32_t rir = hi ? rirLo : rirHi;
422                         if(anchorOverhang > 0) {
423                                 // Does the non-anchor part of the alignment (the
424                                 // "overhang") ruin it?
425                                 bool skipCandidate = false;
426                                 for(uint32_t j = 0; j < anchorOverhang; j++) {
427                                         assert_lt(ri + anchorBitPairs + j, end);
428                                         int rc = (int)ref[rir + anchorBitPairs + j];
429                                         if(rc == 4) {
430                                                 // Oops, encountered an N in the reference in
431                                                 // the overhang portion of the candidate
432                                                 // alignment
433                                                 // (Note that we inverted hi earlier)
434                                                 if(hi) {
435                                                         // Right-to-left
436                                                         skipRightToLefts = anchorOverhang - j - 1;
437                                                 } else {
438                                                         // Left-to-right
439                                                         skipLeftToRights = anchorBitPairs + j;
440                                                 }
441                                                 skipCandidate = true;
442                                                 break; // Skip this candidate
443                                         }
444                                         if((int)qry[32 + j] != rc) {
445                                                 // Yes, overhang ruins it
446                                                 foundHit = false;
447                                                 break;
448                                         }
449                                 }
450                                 if(skipCandidate) continue;
451                         }
452                         if(foundHit) {
453                                 if(pairs != NULL) {
454                                         TU64Pair p;
455                                         if(ri < aoff) {
456                                                 // By convention, the upstream mate's
457                                                 // coordinates go in the 'first' field
458                                                 p.first  = ((uint64_t)tidx << 32) | (uint64_t)ri;
459                                                 p.second = ((uint64_t)tidx << 32) | (uint64_t)aoff;
460                                         } else {
461                                                 p.first  = ((uint64_t)tidx << 32) | (uint64_t)aoff;
462                                                 p.second = ((uint64_t)tidx << 32) | (uint64_t)ri;
463                                         }
464                                         if(pairs->find(p) != pairs->end()) {
465                                                 // We already found this hit!  Continue.
466                                                 ASSERT_ONLY(duplicates++);
467                                                 ASSERT_ONLY(r2i++);
468                                                 continue;
469                                         } else {
470                                                 // Record this hit
471                                                 pairs->insert(p);
472                                         }
473                                 }
474                                 assert_lt(r2i, r2.size());
475                                 assert_eq(re2[r2i], ri);
476                                 ranges.resize(ranges.size()+1);
477                                 Range& range = ranges.back();
478                                 range.stratum = 0;
479                                 range.numMms = 0;
480                                 assert_eq(0, range.mms.size());
481                                 assert_eq(0, range.refcs.size());
482                                 ASSERT_ONLY(r2i++);
483                                 results.push_back(ri);
484                                 if(--numToFind == 0) return;
485                         } else {
486                                 // Keep scanning
487                         }
488                 }
489                 assert_leq(duplicates, r2.size());
490                 assert_geq(r2.size() - duplicates, ranges.size() - rangesInitSz);
491                 return; // no hit
492         }
493 };
494
495 /**
496  * Defined in ref_aligner.cpp.  Maps an octet representing the XOR of
497  * two two-bit-per-base-encoded DNA sequences to the number of bases
498  * that mismatch between the two.
499  */
500 extern unsigned char u8toMms[];
501
502 /**
503  * Concrete RefAligner for finding nearby 1-mismatch hits given an
504  * anchor hit.
505  */
506 template<typename TStr>
507 class OneMMRefAligner : public RefAligner<TStr> {
508
509         typedef seqan::String<seqan::Dna5> TDna5Str;
510         typedef seqan::String<char> TCharStr;
511         typedef std::vector<Range> TRangeVec;
512         typedef std::vector<uint32_t> TU32Vec;
513         typedef std::pair<uint64_t, uint64_t> TU64Pair;
514         typedef std::set<TU64Pair> TSetPairs;
515
516 public:
517
518         OneMMRefAligner(bool color, bool verbose, bool quiet) :
519                 RefAligner<TStr>(color, verbose, quiet) { }
520
521         virtual ~OneMMRefAligner() { }
522
523 protected:
524         /**
525          * Because we're doing end-to-end exact, we don't care which end of
526          * 'qry' is the 5' end.
527          */
528         void naiveFind(uint32_t numToFind,
529                        uint32_t tidx,
530                        uint8_t* ref,
531                    const TDna5Str& qry,
532                    const TCharStr& quals,
533                    uint32_t begin,
534                    uint32_t end,
535                        TRangeVec& ranges,
536                        TU32Vec& results,
537                        TSetPairs* pairs,
538                        uint32_t aoff,
539                        bool seedOnLeft) const
540         {
541                 assert_gt(numToFind, 0);
542                 const uint32_t qlen = seqan::length(qry);
543                 assert_geq(end - begin, qlen); // caller should have checked this
544                 assert_gt(end, begin);
545                 assert_gt(qlen, 0);
546                 uint32_t qend = end - qlen;
547                 uint32_t lim = qend - begin;
548                 uint32_t halfway = begin + (lim >> 1);
549                 bool hi = false;
550                 for(uint32_t i = 1; i <= lim+1; i++) {
551                         uint32_t ri;  // leftmost position in candidate alignment
552                         uint32_t rir; // same, minus begin; for indexing into ref[]
553                         if(hi) {
554                                 ri = halfway + (i >> 1); rir = ri - begin;
555                                 assert_leq(ri, qend);
556                         } else {
557                                 ri = halfway - (i >> 1); rir = ri - begin;
558                                 assert_geq(ri, begin);
559                         }
560                         hi = !hi;
561                         // Do the naive comparison
562                         bool match = true;
563                         int refc = -1;
564                         uint32_t mmOff = 0xffffffff;
565                         int mms = 0;
566                         for(uint32_t j = 0; j < qlen; j++) {
567 #if 0
568                                 // Count Ns in the reference as mismatches
569                                 const int q = (int)qry[j];
570                                 const int r = (int)ref[rir + j];
571                                 assert_leq(q, 4);
572                                 assert_leq(r, 4);
573                                 if(q == 4 || r == 4 || q != r) {
574 #else
575                                 // Disallow alignments that involve an N in the
576                                 // reference
577                                 const int r = (int)ref[rir + j];
578                                 if(r & 4) {
579                                         // N in reference; bail
580                                         match = false;
581                                         break;
582                                 }
583                                 const int q = (int)qry[j];
584                                 assert_leq(q, 4);
585                                 assert_lt(r, 4);
586                                 if(q != r) {
587 #endif
588                                         // Mismatch!
589                                         if(++mms > 1) {
590                                                 // Too many; reject this alignment
591                                                 match = false;
592                                                 break;
593                                         } else {
594                                                 // First one; remember offset and ref char
595                                                 refc = "ACGTN"[r];
596                                                 mmOff = j;
597                                         }
598                                 }
599                         }
600                         if(match) {
601                                 assert_leq(mms, 1);
602                                 ranges.resize(ranges.size()+1);
603                                 Range& range = ranges.back();
604                                 range.stratum = mms;
605                                 range.numMms = mms;
606                                 assert_eq(0, range.mms.size());
607                                 assert_eq(0, range.refcs.size());
608                                 if(mms == 1) {
609                                         assert_lt(mmOff, qlen);
610                                         range.mms.push_back(mmOff);
611                                         range.refcs.push_back(refc);
612                                 }
613                                 results.push_back(ri);
614                         }
615                 }
616                 return;
617         }
618
619         /**
620          * Because we're doing end-to-end exact, we don't care which end of
621          * 'qry' is the 5' end.
622          */
623         virtual void anchor64Find(uint32_t numToFind,
624                         uint32_t tidx,
625                         uint8_t* ref,
626                     const TDna5Str& qry,
627                     const TCharStr& quals,
628                     uint32_t begin,
629                     uint32_t end,
630                         TRangeVec& ranges,
631                         TU32Vec& results,
632                         TSetPairs* pairs = NULL,
633                         uint32_t aoff = 0xffffffff,
634                         bool seedOnLeft = false) const
635         {
636                 assert_gt(numToFind, 0);
637                 ASSERT_ONLY(const uint32_t rangesInitSz = ranges.size());
638                 ASSERT_ONLY(uint32_t duplicates = 0);
639                 ASSERT_ONLY(uint32_t r2i = 0);
640                 const uint32_t qlen = seqan::length(qry);
641                 assert_geq(end - begin, qlen); // caller should have checked this
642                 assert_gt(end, begin);
643                 assert_gt(qlen, 0);
644 #ifndef NDEBUG
645                 // Get results from the naive matcher for sanity-checking
646                 TRangeVec r2; TU32Vec re2;
647                 naiveFind(numToFind, tidx, ref, qry, quals, begin, end, r2,
648                           re2, pairs, aoff, seedOnLeft);
649 #endif
650                 const uint32_t anchorBitPairs = min<int>(qlen, 32);
651                 const int lhsShift = ((anchorBitPairs - 1) << 1);
652                 const uint32_t anchorCushion  = 32 - anchorBitPairs;
653                 // anchorOverhang = # read bases not included in the anchor
654                 const uint32_t anchorOverhang = (qlen <= 32 ? 0 : (qlen - 32));
655                 const uint32_t lim = end - qlen - begin;
656                 const uint32_t halfway = begin + (lim >> 1);
657                 uint64_t anchor = 0llu;
658                 uint64_t buffw = 0llu;  // rotating ref sequence buffer
659                 // OR the 'diff' buffer with this so that we can always count
660                 // 'N's as mismatches
661                 uint64_t diffMask = 0llu;
662                 // Set up a mask that we'll apply to the two bufs every round
663                 // to discard bits that were rotated out of the anchor area
664                 uint64_t clearMask = 0xffffffffffffffffllu;
665                 bool useMask = false;
666                 if(anchorBitPairs < 32) {
667                         clearMask >>= ((32-anchorBitPairs) << 1);
668                         useMask = true;
669                 }
670                 int nsInAnchor = 0;
671                 int nPos = -1;
672                 uint32_t skipLeftToRights = 0;
673                 uint32_t skipRightToLefts = 0;
674                 // Construct the 'anchor' 64-bit buffer so that it holds all of
675                 // the first 'anchorBitPairs' bit pairs of the query.
676                 for(uint32_t i = 0; i < anchorBitPairs; i++) {
677                         int c = (int)qry[i]; // next query character
678                         int r = (int)ref[halfway - begin + i]; // next reference character
679                         if(r & 4) {
680                                 r = 0;
681                                 // The right-to-left direction absorbs the candidate
682                                 // alignment based at 'halfway'; so skipLeftToRights is
683                                 // i, not i+1
684                                 skipLeftToRights = max(skipLeftToRights, i);
685                                 skipRightToLefts = max(skipRightToLefts, anchorBitPairs - i);
686                         }
687                         assert_leq(c, 4);
688                         assert_lt(r, 4);
689                         // Special case: query has an 'N'
690                         if(c == 4) {
691                                 if(++nsInAnchor > 1) {
692                                         // More than one 'N' in the anchor region; can't
693                                         // possibly have a 1-mismatch hit anywhere
694                                         assert_eq(r2.size(), ranges.size() - rangesInitSz);
695                                         return;   // can't match if query has Ns
696                                 }
697                                 nPos = (int)i;
698                                 // Make it look like an 'A' in the anchor
699                                 c = 0;
700                                 diffMask = (diffMask << 2llu) | 1llu;
701                         } else {
702                                 diffMask <<= 2llu;
703                         }
704                         anchor  = ((anchor  << 2llu) | c);
705                         buffw = ((buffw << 2llu) | r);
706                 }
707                 // Check whether read is disqualified by Ns outside of the anchor
708                 // region
709                 for(uint32_t i = anchorBitPairs; i < qlen; i++) {
710                         if((int)qry[i] == 4) {
711                                 if(++nsInAnchor > 1) {
712                                         assert_eq(r2.size(), ranges.size() - rangesInitSz);
713                                         return; // can't match if query has Ns
714                                 }
715                         }
716                 }
717                 uint64_t bufbw = buffw;
718                 // We're moving the right-hand edge of the anchor along until
719                 // it's 'anchorOverhang' chars from the end of the target region.
720                 // Note that we're not making a 3'/5' distinction here; if we
721                 // were, we might need to make the 'anchorOverhang' adjustment on
722                 // the left end of the range rather than the right.
723                 bool hi = false;
724                 uint32_t riHi  = halfway;
725                 uint32_t rirHi = halfway - begin;
726                 uint32_t rirHiAnchor = rirHi + anchorBitPairs - 1;
727                 uint32_t riLo  = halfway + 1;
728                 uint32_t rirLo = halfway - begin + 1;
729                 for(uint32_t i = 1; i <= lim + 1; i++) {
730                         int r;       // new reference char
731                         uint64_t diff;
732                         assert_lt(skipLeftToRights, qlen);
733                         assert_leq(skipRightToLefts, qlen);
734                         if(hi) {
735                                 hi = false;
736                                 // Moving left-to-right
737                                 riHi++;
738                                 rirHi++;
739                                 rirHiAnchor++;
740                                 r = (int)ref[rirHiAnchor];
741                                 if(r & 4) {
742                                         r = 0;
743                                         skipLeftToRights = anchorBitPairs;
744                                 }
745                                 assert_lt(r, 4);
746                                 // Bring in new base pair at the least significant
747                                 // position
748                                 buffw = ((buffw << 2llu) | r);
749                                 if(useMask) buffw &= clearMask;
750                                 if(skipLeftToRights > 0) {
751                                         skipLeftToRights--;
752                                         continue;
753                                 }
754                                 diff = (buffw ^ anchor) | diffMask;
755                         } else {
756                                 hi = true;
757                                 // Moving right-to-left
758                                 riLo--;
759                                 rirLo--;
760                                 r = (int)ref[rirLo];
761                                 if(r & 4) {
762                                         r = 0;
763                                         skipRightToLefts = qlen;
764                                 }
765                                 assert_lt(r, 4);
766                                 if(i >= 2) {
767                                         bufbw >>= 2llu;
768                                         // Bring in new base pair at the most significant
769                                         // position
770                                         bufbw |= ((uint64_t)r << lhsShift);
771                                 }
772                                 if(skipRightToLefts > 0) {
773                                         skipRightToLefts--;
774                                         continue;
775                                 }
776                                 diff = (bufbw ^ anchor) | diffMask;
777                         }
778                         if((diff & 0xffffffff00000000llu) &&
779                            (diff & 0x00000000ffffffffllu)) continue;
780                         uint32_t ri  = hi ? riLo  : riHi;
781                         uint32_t rir = hi ? rirLo : rirHi;
782                         // Could use pop count
783                         uint8_t *diff8 = reinterpret_cast<uint8_t*>(&diff);
784                         // As a first cut, see if there are too many mismatches in
785                         // the first and last parts of the anchor
786                         uint32_t diffs = u8toMms[(int)diff8[0]] + u8toMms[(int)diff8[7]];
787                         if(diffs > 1) continue;
788                         diffs += u8toMms[(int)diff8[1]] +
789                                  u8toMms[(int)diff8[2]] +
790                                  u8toMms[(int)diff8[3]] +
791                                  u8toMms[(int)diff8[4]] +
792                                  u8toMms[(int)diff8[5]] +
793                                  u8toMms[(int)diff8[6]];
794                         uint32_t mmpos = 0xffffffff;
795                         int refc = -1;
796                         if(diffs > 1) {
797                                 // Too many differences
798                                 continue;
799                         } else if(diffs == 1 && nPos != -1) {
800                                 // There was one difference, but there was also one N,
801                                 // so we already know where the difference is
802                                 mmpos = nPos;
803                                 refc = "ACGT"[(int)ref[rir + nPos]];
804                         } else if(diffs == 1) {
805                                 // Figure out which position mismatched
806                                 mmpos = 31;
807                                 if((diff & 0xffffffffllu) == 0) { diff >>= 32llu; mmpos -= 16; }
808                                 assert_neq(0, diff);
809                                 if((diff & 0xffffllu) == 0)     { diff >>= 16llu; mmpos -=  8; }
810                                 assert_neq(0, diff);
811                                 if((diff & 0xffllu) == 0)       { diff >>= 8llu;  mmpos -=  4; }
812                                 assert_neq(0, diff);
813                                 if((diff & 0xfllu) == 0)        { diff >>= 4llu;  mmpos -=  2; }
814                                 assert_neq(0, diff);
815                                 if((diff & 0x3llu) == 0)        { mmpos--; }
816                                 assert_neq(0, diff);
817                                 assert_geq(mmpos, 0);
818                                 assert_lt(mmpos, 32);
819                                 mmpos -= anchorCushion;
820                                 refc = "ACGT"[(int)ref[rir + mmpos]];
821                         }
822                         // Now extend the anchor into a longer alignment
823                         bool foundHit = true;
824                         if(anchorOverhang > 0) {
825                                 assert_leq(ri + anchorBitPairs + anchorOverhang, end);
826                                 bool skipCandidate = false;
827                                 for(uint32_t j = 0; j < anchorOverhang; j++) {
828                                         int rc = (int)ref[rir + 32 + j];
829                                         if(rc == 4) {
830                                                 // Oops, encountered an N in the reference in
831                                                 // the overhang portion of the candidate
832                                                 // alignment
833                                                 // (Note that we inverted hi earlier)
834                                                 if(hi) {
835                                                         // Right-to-left
836                                                         skipRightToLefts = anchorOverhang - j - 1;
837                                                 } else {
838                                                         // Left-to-right
839                                                         skipLeftToRights = anchorBitPairs + j;
840                                                 }
841                                                 skipCandidate = true; // Skip this candidate
842                                                 break;
843                                         }
844                                         if((int)qry[32 + j] != rc) {
845                                                 if(++diffs > 1) {
846                                                         foundHit = false;
847                                                         break;
848                                                 } else {
849                                                         mmpos = 32 + j;
850                                                         refc = "ACGT"[(int)ref[rir + 32 + j]];
851                                                 }
852                                         }
853                                 }
854                                 if(skipCandidate) continue;
855                         }
856                         if(!foundHit) continue;
857                         if(pairs != NULL) {
858                                 TU64Pair p;
859                                 if(ri < aoff) {
860                                         // By convention, the upstream mate's
861                                         // coordinates go in the 'first' field
862                                         p.first  = ((uint64_t)tidx << 32) | (uint64_t)ri;
863                                         p.second = ((uint64_t)tidx << 32) | (uint64_t)aoff;
864                                 } else {
865                                         p.second = ((uint64_t)tidx << 32) | (uint64_t)ri;
866                                         p.first  = ((uint64_t)tidx << 32) | (uint64_t)aoff;
867                                 }
868                                 if(pairs->find(p) != pairs->end()) {
869                                         // We already found this hit!  Continue.
870                                         ASSERT_ONLY(duplicates++);
871                                         ASSERT_ONLY(r2i++);
872                                         continue;
873                                 } else {
874                                         // Record this hit
875                                         pairs->insert(p);
876                                 }
877                         }
878                         assert_leq(diffs, 1);
879                         assert_lt(r2i, r2.size());
880                         assert_eq(re2[r2i], ri);
881                         ranges.resize(ranges.size()+1);
882                         Range& range = ranges.back();
883                         assert_eq(diffs, r2[r2i].numMms);
884                         range.stratum = diffs;
885                         range.numMms = diffs;
886                         assert_eq(0, range.mms.size());
887                         assert_eq(0, range.refcs.size());
888                         if(diffs == 1) {
889                                 assert_neq(mmpos, 0xffffffff);
890                                 assert_eq(mmpos, r2[r2i].mms[0]);
891                                 assert_neq(-1, refc);
892                                 assert_eq(refc, r2[r2i].refcs[0]);
893                                 range.mms.push_back(mmpos);
894                                 range.refcs.push_back(refc);
895                         }
896                         ASSERT_ONLY(r2i++);
897                         results.push_back(ri);
898                         if(--numToFind == 0) return;
899                 }
900                 assert_leq(duplicates, r2.size());
901                 assert_geq(r2.size() - duplicates, ranges.size() - rangesInitSz);
902                 return; // no hit
903         }
904 };
905
906 /**
907  * Concrete RefAligner for finding nearby 2-mismatch hits given an
908  * anchor hit.
909  */
910 template<typename TStr>
911 class TwoMMRefAligner : public RefAligner<TStr> {
912
913         typedef seqan::String<seqan::Dna5> TDna5Str;
914         typedef seqan::String<char> TCharStr;
915         typedef std::vector<Range> TRangeVec;
916         typedef std::vector<uint32_t> TU32Vec;
917         typedef std::pair<uint64_t, uint64_t> TU64Pair;
918         typedef std::set<TU64Pair> TSetPairs;
919
920 public:
921
922         TwoMMRefAligner(bool color, bool verbose, bool quiet) :
923                 RefAligner<TStr>(color, verbose, quiet) { }
924
925         virtual ~TwoMMRefAligner() { }
926
927 protected:
928         /**
929          * Because we're doing end-to-end exact, we don't care which end of
930          * 'qry' is the 5' end.
931          */
932         void naiveFind(uint32_t numToFind,
933                                    uint32_t tidx,
934                                    uint8_t* ref,
935                                    const TDna5Str& qry,
936                                    const TCharStr& quals,
937                                    uint32_t begin,
938                                    uint32_t end,
939                                    TRangeVec& ranges,
940                                    TU32Vec& results,
941                                    TSetPairs* pairs,
942                                    uint32_t aoff,
943                                    bool seedOnLeft) const
944         {
945                 assert_gt(numToFind, 0);
946                 const uint32_t qlen = seqan::length(qry);
947                 assert_geq(end - begin, qlen); // caller should have checked this
948                 assert_gt(end, begin);
949                 assert_gt(qlen, 0);
950                 uint32_t qend = end - qlen;
951                 uint32_t lim = qend - begin;
952                 uint32_t halfway = begin + (lim >> 1);
953                 bool hi = false;
954                 for(uint32_t i = 1; i <= lim+1; i++) {
955                         uint32_t ri;  // leftmost position in candidate alignment
956                         uint32_t rir; // same, minus begin; for indexing into ref[]
957                         if(hi) {
958                                 ri = halfway + (i >> 1); rir = ri - begin;
959                                 assert_leq(ri, qend);
960                         } else {
961                                 ri = halfway - (i >> 1); rir = ri - begin;
962                                 assert_geq(ri, begin);
963                         }
964                         hi = !hi;
965                         // Do the naive comparison
966                         bool match = true;
967                         int refc1 = -1;
968                         uint32_t mmOff1 = 0xffffffff;
969                         int refc2 = -1;
970                         uint32_t mmOff2 = 0xffffffff;
971                         int mms = 0;
972                         for(uint32_t j = 0; j < qlen; j++) {
973 #if 0
974                                 // Count Ns in the reference as mismatches
975                                 const int q = (int)qry[j];
976                                 const int r = (int)ref[rir + j];
977                                 assert_leq(q, 4);
978                                 assert_leq(r, 4);
979                                 if(q == 4 || r == 4 || q != r) {
980 #else
981                                 // Disallow alignments that involve an N in the
982                                 // reference
983                                 const int r = (int)ref[rir + j];
984                                 if(r & 4) {
985                                         // N in reference; bail
986                                         match = false;
987                                         break;
988                                 }
989                                 const int q = (int)qry[j];
990                                 assert_leq(q, 4);
991                                 assert_lt(r, 4);
992                                 if(q != r) {
993 #endif
994                                         // Mismatch!
995                                         if(++mms > 2) {
996                                                 // Too many; reject this alignment
997                                                 match = false;
998                                                 break;
999                                         } else if(mms == 2) {
1000                                                 // Second one; remember offset and ref char
1001                                                 refc2 = "ACGTN"[r];
1002                                                 mmOff2 = j;
1003                                         } else {
1004                                                 assert_eq(1, mms);
1005                                                 // First one; remember offset and ref char
1006                                                 refc1 = "ACGTN"[r];
1007                                                 mmOff1 = j;
1008                                         }
1009                                 }
1010                         }
1011                         if(match) {
1012                                 assert_leq(mms, 2);
1013                                 ranges.resize(ranges.size()+1);
1014                                 Range& range = ranges.back();
1015                                 range.stratum = mms;
1016                                 range.numMms = mms;
1017                                 assert_eq(0, range.mms.size());
1018                                 assert_eq(0, range.refcs.size());
1019                                 if(mms > 0) {
1020                                         assert_lt(mmOff1, qlen);
1021                                         assert(refc1 == 'A' || refc1 == 'C' || refc1 == 'G' || refc1 == 'T');
1022                                         range.mms.push_back(mmOff1);
1023                                         range.refcs.push_back(refc1);
1024                                         if(mms > 1) {
1025                                                 assert_eq(2, mms);
1026                                                 assert_lt(mmOff2, qlen);
1027                                                 assert(refc2 == 'A' || refc2 == 'C' || refc2 == 'G' || refc2 == 'T');
1028                                                 range.mms.push_back(mmOff2);
1029                                                 range.refcs.push_back(refc2);
1030                                         }
1031                                 }
1032                                 results.push_back(ri);
1033                         }
1034                 }
1035                 return;
1036         }
1037
1038         /**
1039          * Because we're doing end-to-end exact, we don't care which end of
1040          * 'qry' is the 5' end.
1041          */
1042         virtual void anchor64Find(uint32_t numToFind,
1043                                         uint32_t tidx,
1044                                         uint8_t* ref,
1045                                         const TDna5Str& qry,
1046                                         const TCharStr& quals,
1047                                         uint32_t begin,
1048                                         uint32_t end,
1049                                         TRangeVec& ranges,
1050                                         TU32Vec& results,
1051                                         TSetPairs* pairs = NULL,
1052                                         uint32_t aoff = 0xffffffff,
1053                                         bool seedOnLeft = false) const
1054         {
1055                 assert_gt(numToFind, 0);
1056                 ASSERT_ONLY(const uint32_t rangesInitSz = ranges.size());
1057                 ASSERT_ONLY(uint32_t duplicates = 0);
1058                 ASSERT_ONLY(uint32_t r2i = 0);
1059                 const uint32_t qlen = seqan::length(qry);
1060                 assert_geq(end - begin, qlen); // caller should have checked this
1061                 assert_gt(end, begin);
1062                 assert_gt(qlen, 0);
1063 #ifndef NDEBUG
1064                 // Get results from the naive matcher for sanity-checking
1065                 TRangeVec r2; TU32Vec re2;
1066                 naiveFind(numToFind, tidx, ref, qry, quals, begin, end, r2,
1067                                   re2, pairs, aoff, seedOnLeft);
1068 #endif
1069                 const uint32_t anchorBitPairs = min<int>(qlen, 32);
1070                 const int lhsShift = ((anchorBitPairs - 1) << 1);
1071                 const uint32_t anchorCushion  = 32 - anchorBitPairs;
1072                 // anchorOverhang = # read bases not included in the anchor
1073                 const uint32_t anchorOverhang = (qlen <= 32 ? 0 : (qlen - 32));
1074                 const uint32_t lim = end - qlen - begin;
1075                 const uint32_t halfway = begin + (lim >> 1);
1076                 uint64_t anchor = 0llu;
1077                 uint64_t buffw = 0llu;  // rotating ref sequence buffer
1078                 // OR the 'diff' buffer with this so that we can always count
1079                 // 'N's as mismatches
1080                 uint64_t diffMask = 0llu;
1081                 // Set up a mask that we'll apply to the two bufs every round
1082                 // to discard bits that were rotated out of the anchor area
1083                 uint64_t clearMask = 0xffffffffffffffffllu;
1084                 bool useMask = false;
1085                 if(anchorBitPairs < 32) {
1086                         clearMask >>= ((32-anchorBitPairs) << 1);
1087                         useMask = true;
1088                 }
1089                 int nsInAnchor = 0;
1090                 uint32_t nPoss = 0;
1091                 int nPos1 = -1;
1092                 int nPos2 = -1;
1093                 uint32_t skipLeftToRights = 0;
1094                 uint32_t skipRightToLefts = 0;
1095                 // Construct the 'anchor' 64-bit buffer so that it holds all of
1096                 // the first 'anchorBitPairs' bit pairs of the query.
1097                 for(uint32_t i = 0; i < anchorBitPairs; i++) {
1098                         int c = (int)qry[i]; // next query character
1099                         int r = (int)ref[halfway - begin + i]; // next reference character
1100                         if(r & 4) {
1101                                 r = 0;
1102                                 // The right-to-left direction absorbs the candidate
1103                                 // alignment based at 'halfway'; so skipLeftToRights is
1104                                 // i, not i+1
1105                                 skipLeftToRights = max(skipLeftToRights, i);
1106                                 skipRightToLefts = max(skipRightToLefts, anchorBitPairs - i);
1107                         }
1108                         assert_leq(c, 4);
1109                         assert_lt(r, 4);
1110                         // Special case: query has an 'N'
1111                         if(c == 4) {
1112                                 if(++nsInAnchor > 2) {
1113                                         // More than two 'N's in the anchor region; can't
1114                                         // possibly have a 2-mismatch hit anywhere
1115                                         return;   // can't match if query has Ns
1116                                 } else if(nsInAnchor == 2) {
1117                                         nPos2 = (int)i;
1118                                         nPoss++;
1119                                 } else {
1120                                         assert_eq(1, nsInAnchor);
1121                                         nPos1 = (int)i;
1122                                         nPoss++;
1123                                 }
1124                                 // Make it look like an 'A' in the anchor
1125                                 c = 0;
1126                                 diffMask = (diffMask << 2llu) | 1llu;
1127                         } else {
1128                                 diffMask <<= 2llu;
1129                         }
1130                         anchor  = ((anchor  << 2llu) | c);
1131                         buffw = ((buffw << 2llu) | r);
1132                 }
1133                 assert_leq(nPoss, 2);
1134                 // Check whether read is disqualified by Ns outside of the anchor
1135                 // region
1136                 for(uint32_t i = anchorBitPairs; i < qlen; i++) {
1137                         if((int)qry[i] == 4) {
1138                                 if(++nsInAnchor > 2) {
1139                                         return; // can't match if query has Ns
1140                                 }
1141                         }
1142                 }
1143                 uint64_t bufbw = buffw;
1144                 // We're moving the right-hand edge of the anchor along until
1145                 // it's 'anchorOverhang' chars from the end of the target region.
1146                 // Note that we're not making a 3'/5' distinction here; if we
1147                 // were, we might need to make the 'anchorOverhang' adjustment on
1148                 // the left end of the range rather than the right.
1149                 bool hi = false;
1150                 uint32_t riHi  = halfway;
1151                 uint32_t rirHi = halfway - begin;
1152                 uint32_t rirHiAnchor = rirHi + anchorBitPairs - 1;
1153                 uint32_t riLo  = halfway + 1;
1154                 uint32_t rirLo = halfway - begin + 1;
1155                 uint32_t i;
1156                 for(i = 1; i <= lim + 1; i++) {
1157                         int r;       // new reference char
1158                         uint64_t diff;
1159                         assert_lt(skipLeftToRights, qlen);
1160                         assert_leq(skipRightToLefts, qlen);
1161                         if(hi) {
1162                                 hi = false;
1163                                 // Moving left-to-right
1164                                 riHi++;
1165                                 rirHi++;
1166                                 rirHiAnchor++;
1167                                 r = (int)ref[rirHiAnchor];
1168                                 if(r & 4) {
1169                                         r = 0;
1170                                         skipLeftToRights = anchorBitPairs;
1171                                 }
1172                                 assert_lt(r, 4);
1173                                 // Bring in new base pair at the least significant
1174                                 // position
1175                                 buffw = ((buffw << 2llu) | r);
1176                                 if(useMask) buffw &= clearMask;
1177                                 if(skipLeftToRights > 0) {
1178                                         skipLeftToRights--;
1179                                         continue;
1180                                 }
1181                                 diff = (buffw ^ anchor) | diffMask;
1182                         } else {
1183                                 hi = true;
1184                                 // Moving right-to-left
1185                                 riLo--;
1186                                 rirLo--;
1187                                 r = (int)ref[rirLo];
1188                                 if(r & 4) {
1189                                         r = 0;
1190                                         skipRightToLefts = qlen;
1191                                 }
1192                                 assert_lt(r, 4);
1193                                 if(i >= 2) {
1194                                         bufbw >>= 2llu;
1195                                         // Bring in new base pair at the most significant
1196                                         // position
1197                                         bufbw |= ((uint64_t)r << lhsShift);
1198                                 }
1199                                 if(skipRightToLefts > 0) {
1200                                         skipRightToLefts--;
1201                                         continue;
1202                                 }
1203                                 diff = (bufbw ^ anchor) | diffMask;
1204                         }
1205                         if((diff & 0xfffff00000000000llu) &&
1206                            (diff & 0x00000ffffff00000llu) &&
1207                            (diff & 0x00000000000fffffllu)) continue;
1208                         uint32_t ri  = hi ? riLo  : riHi;
1209                         uint32_t rir = hi ? rirLo : rirHi;
1210                         // Could use pop count
1211                         uint8_t *diff8 = reinterpret_cast<uint8_t*>(&diff);
1212                         // As a first cut, see if there are too many mismatches in
1213                         // the first and last parts of the anchor
1214                         uint32_t diffs = u8toMms[(int)diff8[0]] + u8toMms[(int)diff8[7]];
1215                         if(diffs > 2) continue;
1216                         diffs += u8toMms[(int)diff8[1]] +
1217                                          u8toMms[(int)diff8[2]] +
1218                                          u8toMms[(int)diff8[3]] +
1219                                          u8toMms[(int)diff8[4]] +
1220                                          u8toMms[(int)diff8[5]] +
1221                                          u8toMms[(int)diff8[6]];
1222                         uint32_t mmpos1 = 0xffffffff;
1223                         int refc1 = -1;
1224                         uint32_t mmpos2 = 0xffffffff;
1225                         int refc2 = -1;
1226                         if(diffs > 2) {
1227                                 // Too many differences
1228                                 continue;
1229                         } else if(nPoss > 1 && diffs == nPoss) {
1230                                 // There was one difference, but there was also one N,
1231                                 // so we already know where the difference is
1232                                 mmpos1 = nPos1;
1233                                 refc1 = "ACGT"[(int)ref[rir + nPos1]];
1234                                 if(nPoss == 2) {
1235                                         mmpos2 = nPos2;
1236                                         refc2 = "ACGT"[(int)ref[rir + nPos2]];
1237                                 }
1238                         } else if(diffs > 0) {
1239                                 // Figure out which position mismatched
1240                                 uint64_t diff2 = diff;
1241                                 mmpos1 = 31;
1242                                 if((diff & 0xffffffffllu) == 0) { diff >>= 32llu; mmpos1 -= 16; }
1243                                 assert_neq(0, diff);
1244                                 if((diff & 0xffffllu) == 0)     { diff >>= 16llu; mmpos1 -=  8; }
1245                                 assert_neq(0, diff);
1246                                 if((diff & 0xffllu) == 0)       { diff >>= 8llu;  mmpos1 -=  4; }
1247                                 assert_neq(0, diff);
1248                                 if((diff & 0xfllu) == 0)        { diff >>= 4llu;  mmpos1 -=  2; }
1249                                 assert_neq(0, diff);
1250                                 if((diff & 0x3llu) == 0)        { mmpos1--; }
1251                                 assert_neq(0, diff);
1252                                 assert_geq(mmpos1, 0);
1253                                 assert_lt(mmpos1, 32);
1254                                 mmpos1 -= anchorCushion;
1255                                 refc1 = "ACGT"[(int)ref[rir + mmpos1]];
1256                                 if(diffs > 1) {
1257                                         // Figure out the second mismatched position
1258                                         ASSERT_ONLY(uint64_t origDiff2 = diff2);
1259                                         diff2 &= ~(0xc000000000000000llu >> (uint64_t)((mmpos1+anchorCushion) << 1));
1260                                         assert_neq(diff2, origDiff2);
1261                                         mmpos2 = 31;
1262                                         if((diff2 & 0xffffffffllu) == 0) { diff2 >>= 32llu; mmpos2 -= 16; }
1263                                         assert_neq(0, diff2);
1264                                         if((diff2 & 0xffffllu) == 0)     { diff2 >>= 16llu; mmpos2 -=  8; }
1265                                         assert_neq(0, diff2);
1266                                         if((diff2 & 0xffllu) == 0)       { diff2 >>= 8llu;  mmpos2 -=  4; }
1267                                         assert_neq(0, diff2);
1268                                         if((diff2 & 0xfllu) == 0)        { diff2 >>= 4llu;  mmpos2 -=  2; }
1269                                         assert_neq(0, diff2);
1270                                         if((diff2 & 0x3llu) == 0)        { mmpos2--; }
1271                                         assert_neq(0, diff2);
1272                                         assert_geq(mmpos2, 0);
1273                                         assert_lt(mmpos2, 32);
1274                                         mmpos2 -= anchorCushion;
1275                                         assert_neq(mmpos1, mmpos2);
1276                                         refc2 = "ACGT"[(int)ref[rir + mmpos2]];
1277                                         if(mmpos2 < mmpos1) {
1278                                                 uint32_t mmtmp = mmpos1;
1279                                                 mmpos1 = mmpos2;
1280                                                 mmpos2 = mmtmp;
1281                                                 int refctmp = refc1;
1282                                                 refc1 = refc2;
1283                                                 refc2 = refctmp;
1284                                         }
1285                                         assert_lt(mmpos1, mmpos2);
1286                                 }
1287                         }
1288                         // Now extend the anchor into a longer alignment
1289                         bool foundHit = true;
1290                         if(anchorOverhang > 0) {
1291                                 assert_leq(ri + anchorBitPairs + anchorOverhang, end);
1292                                 bool skipCandidate = false;
1293                                 for(uint32_t j = 0; j < anchorOverhang; j++) {
1294                                         int rc = (int)ref[rir + 32 + j];
1295                                         if(rc == 4) {
1296                                                 // Oops, encountered an N in the reference in
1297                                                 // the overhang portion of the candidate
1298                                                 // alignment
1299                                                 // (Note that we inverted hi earlier)
1300                                                 if(hi) {
1301                                                         // Right-to-left
1302                                                         skipRightToLefts = anchorOverhang - j - 1;
1303                                                 } else {
1304                                                         // Left-to-right
1305                                                         skipLeftToRights = anchorBitPairs + j;
1306                                                 }
1307                                                 skipCandidate = true; // Skip this candidate
1308                                                 break;
1309                                         }
1310                                         if((int)qry[32 + j] != rc) {
1311                                                 if(++diffs > 2) {
1312                                                         foundHit = false;
1313                                                         break;
1314                                                 } else if(diffs == 2) {
1315                                                         mmpos2 = 32 + j;
1316                                                         refc2 = "ACGT"[(int)ref[rir + 32 + j]];
1317                                                 } else {
1318                                                         assert_eq(1, diffs);
1319                                                         mmpos1 = 32 + j;
1320                                                         refc1 = "ACGT"[(int)ref[rir + 32 + j]];
1321                                                 }
1322                                         }
1323                                 }
1324                                 if(skipCandidate) continue;
1325                         }
1326                         if(!foundHit) continue;
1327                         if(pairs != NULL) {
1328                                 TU64Pair p;
1329                                 if(ri < aoff) {
1330                                         // By convention, the upstream mate's
1331                                         // coordinates go in the 'first' field
1332                                         p.first  = ((uint64_t)tidx << 32) | (uint64_t)ri;
1333                                         p.second = ((uint64_t)tidx << 32) | (uint64_t)aoff;
1334                                 } else {
1335                                         p.second = ((uint64_t)tidx << 32) | (uint64_t)ri;
1336                                         p.first  = ((uint64_t)tidx << 32) | (uint64_t)aoff;
1337                                 }
1338                                 if(pairs->find(p) != pairs->end()) {
1339                                         // We already found this hit!  Continue.
1340                                         ASSERT_ONLY(duplicates++);
1341                                         ASSERT_ONLY(r2i++);
1342                                         continue;
1343                                 } else {
1344                                         // Record this hit
1345                                         pairs->insert(p);
1346                                 }
1347                         }
1348                         assert_leq(diffs, 2);
1349                         assert_lt(r2i, r2.size());
1350                         assert_eq(re2[r2i], ri);
1351                         ranges.resize(ranges.size()+1);
1352                         Range& range = ranges.back();
1353                         assert_eq(diffs, r2[r2i].numMms);
1354                         range.stratum = diffs;
1355                         range.numMms = diffs;
1356                         assert_eq(0, range.mms.size());
1357                         assert_eq(0, range.refcs.size());
1358                         if(diffs > 0) {
1359                                 assert_neq(mmpos1, 0xffffffff);
1360                                 assert_eq(mmpos1, r2[r2i].mms[0]);
1361                                 assert_neq(-1, refc1);
1362                                 assert_eq(refc1, r2[r2i].refcs[0]);
1363                                 range.mms.push_back(mmpos1);
1364                                 range.refcs.push_back(refc1);
1365                                 if(diffs > 1) {
1366                                         assert_neq(mmpos2, 0xffffffff);
1367                                         assert_eq(mmpos2, r2[r2i].mms[1]);
1368                                         assert_neq(-1, refc2);
1369                                         assert_eq(refc2, r2[r2i].refcs[1]);
1370                                         range.mms.push_back(mmpos2);
1371                                         range.refcs.push_back(refc2);
1372                                 }
1373                         }
1374                         ASSERT_ONLY(r2i++);
1375                         results.push_back(ri);
1376                         if(--numToFind == 0) return;
1377                 }
1378                 assert_leq(duplicates, r2.size());
1379                 assert_geq(r2.size() - duplicates, ranges.size() - rangesInitSz);
1380                 return; // no hit
1381         }
1382 };
1383
1384 /**
1385  * Concrete RefAligner for finding nearby 2-mismatch hits given an
1386  * anchor hit.
1387  */
1388 template<typename TStr>
1389 class ThreeMMRefAligner : public RefAligner<TStr> {
1390
1391         typedef seqan::String<seqan::Dna5> TDna5Str;
1392         typedef seqan::String<char> TCharStr;
1393         typedef std::vector<Range> TRangeVec;
1394         typedef std::vector<uint32_t> TU32Vec;
1395         typedef std::pair<uint64_t, uint64_t> TU64Pair;
1396         typedef std::set<TU64Pair> TSetPairs;
1397
1398 public:
1399
1400         ThreeMMRefAligner(bool color, bool verbose, bool quiet) :
1401                 RefAligner<TStr>(color, verbose, quiet) { }
1402
1403         virtual ~ThreeMMRefAligner() { }
1404
1405 protected:
1406         /**
1407          * Because we're doing end-to-end exact, we don't care which end of
1408          * 'qry' is the 5' end.
1409          */
1410         void naiveFind(uint32_t numToFind,
1411                                    uint32_t tidx,
1412                                    uint8_t* ref,
1413                                    const TDna5Str& qry,
1414                                    const TCharStr& quals,
1415                                    uint32_t begin,
1416                                    uint32_t end,
1417                                    TRangeVec& ranges,
1418                                    TU32Vec& results,
1419                                    TSetPairs* pairs,
1420                                    uint32_t aoff,
1421                                    bool seedOnLeft) const
1422         {
1423                 assert_gt(numToFind, 0);
1424                 const uint32_t qlen = seqan::length(qry);
1425                 assert_geq(end - begin, qlen); // caller should have checked this
1426                 assert_gt(end, begin);
1427                 assert_gt(qlen, 0);
1428                 uint32_t qend = end - qlen;
1429                 uint32_t lim = qend - begin;
1430                 uint32_t halfway = begin + (lim >> 1);
1431                 bool hi = false;
1432                 for(uint32_t i = 1; i <= lim+1; i++) {
1433                         uint32_t ri;  // leftmost position in candidate alignment
1434                         uint32_t rir; // same, minus begin; for indexing into ref[]
1435                         if(hi) {
1436                                 ri = halfway + (i >> 1); rir = ri - begin;
1437                                 assert_leq(ri, qend);
1438                         } else {
1439                                 ri = halfway - (i >> 1); rir = ri - begin;
1440                                 assert_geq(ri, begin);
1441                         }
1442                         hi = !hi;
1443                         // Do the naive comparison
1444                         bool match = true;
1445                         int refc1 = -1;
1446                         uint32_t mmOff1 = 0xffffffff;
1447                         int refc2 = -1;
1448                         uint32_t mmOff2 = 0xffffffff;
1449                         int refc3 = -1;
1450                         uint32_t mmOff3 = 0xffffffff;
1451                         int mms = 0;
1452                         for(uint32_t j = 0; j < qlen; j++) {
1453 #if 0
1454                                 // Count Ns in the reference as mismatches
1455                                 const int q = (int)qry[j];
1456                                 const int r = (int)ref[rir + j];
1457                                 assert_leq(q, 4);
1458                                 assert_leq(r, 4);
1459                                 if(q == 4 || r == 4 || q != r) {
1460 #else
1461                                 // Disallow alignments that involve an N in the
1462                                 // reference
1463                                 const int r = (int)ref[rir + j];
1464                                 if(r & 4) {
1465                                         // N in reference; bail
1466                                         match = false;
1467                                         break;
1468                                 }
1469                                 const int q = (int)qry[j];
1470                                 assert_leq(q, 4);
1471                                 assert_lt(r, 4);
1472                                 if(q != r) {
1473 #endif
1474                                         // Mismatch!
1475                                         if(++mms > 3) {
1476                                                 // Too many; reject this alignment
1477                                                 match = false;
1478                                                 break;
1479                                         } else if(mms == 3) {
1480                                                 // Second one; remember offset and ref char
1481                                                 refc3 = "ACGTN"[r];
1482                                                 mmOff3 = j;
1483                                         } else if(mms == 2) {
1484                                                 // Second one; remember offset and ref char
1485                                                 refc2 = "ACGTN"[r];
1486                                                 mmOff2 = j;
1487                                         } else {
1488                                                 assert_eq(1, mms);
1489                                                 // First one; remember offset and ref char
1490                                                 refc1 = "ACGTN"[r];
1491                                                 mmOff1 = j;
1492                                         }
1493                                 }
1494                         }
1495                         if(match) {
1496                                 assert_leq(mms, 3);
1497                                 ranges.resize(ranges.size()+1);
1498                                 Range& range = ranges.back();
1499                                 range.stratum = mms;
1500                                 range.numMms = mms;
1501                                 assert_eq(0, range.mms.size());
1502                                 assert_eq(0, range.refcs.size());
1503                                 if(mms > 0) {
1504                                         assert_lt(mmOff1, qlen);
1505                                         assert(refc1 == 'A' || refc1 == 'C' || refc1 == 'G' || refc1 == 'T');
1506                                         range.mms.push_back(mmOff1);
1507                                         range.refcs.push_back(refc1);
1508                                         if(mms > 1) {
1509                                                 assert_lt(mmOff2, qlen);
1510                                                 assert(refc2 == 'A' || refc2 == 'C' || refc2 == 'G' || refc2 == 'T');
1511                                                 range.mms.push_back(mmOff2);
1512                                                 range.refcs.push_back(refc2);
1513                                                 if(mms > 2) {
1514                                                         assert_eq(3, mms);
1515                                                         assert_lt(mmOff3, qlen);
1516                                                         assert(refc3 == 'A' || refc3 == 'C' || refc3 == 'G' || refc3 == 'T');
1517                                                         range.mms.push_back(mmOff3);
1518                                                         range.refcs.push_back(refc3);
1519                                                 }
1520                                         }
1521                                 }
1522                                 results.push_back(ri);
1523                         }
1524                 }
1525                 return;
1526         }
1527
1528         /**
1529          * Because we're doing end-to-end exact, we don't care which end of
1530          * 'qry' is the 5' end.
1531          */
1532         virtual void anchor64Find(uint32_t numToFind,
1533                                         uint32_t tidx,
1534                                         uint8_t* ref,
1535                                         const TDna5Str& qry,
1536                                         const TCharStr& quals,
1537                                         uint32_t begin,
1538                                         uint32_t end,
1539                                         TRangeVec& ranges,
1540                                         TU32Vec& results,
1541                                         TSetPairs* pairs = NULL,
1542                                         uint32_t aoff = 0xffffffff,
1543                                         bool seedOnLeft = false) const
1544         {
1545                 assert_gt(numToFind, 0);
1546                 ASSERT_ONLY(const uint32_t rangesInitSz = ranges.size());
1547                 ASSERT_ONLY(uint32_t duplicates = 0);
1548                 ASSERT_ONLY(uint32_t r2i = 0);
1549                 const uint32_t qlen = seqan::length(qry);
1550                 assert_geq(end - begin, qlen); // caller should have checked this
1551                 assert_gt(end, begin);
1552                 assert_gt(qlen, 0);
1553 #ifndef NDEBUG
1554                 // Get results from the naive matcher for sanity-checking
1555                 TRangeVec r2; TU32Vec re2;
1556                 naiveFind(numToFind, tidx, ref, qry, quals, begin, end, r2,
1557                                   re2, pairs, aoff, seedOnLeft);
1558 #endif
1559                 const uint32_t anchorBitPairs = min<int>(qlen, 32);
1560                 const int lhsShift = ((anchorBitPairs - 1) << 1);
1561                 const uint32_t anchorCushion  = 32 - anchorBitPairs;
1562                 // anchorOverhang = # read bases not included in the anchor
1563                 const uint32_t anchorOverhang = (qlen <= 32 ? 0 : (qlen - 32));
1564                 const uint32_t lim = end - qlen - begin;
1565                 const uint32_t halfway = begin + (lim >> 1);
1566                 uint64_t anchor = 0llu;
1567                 uint64_t buffw = 0llu;  // rotating ref sequence buffer
1568                 // OR the 'diff' buffer with this so that we can always count
1569                 // 'N's as mismatches
1570                 uint64_t diffMask = 0llu;
1571                 // Set up a mask that we'll apply to the two bufs every round
1572                 // to discard bits that were rotated out of the anchor area
1573                 uint64_t clearMask = 0xffffffffffffffffllu;
1574                 bool useMask = false;
1575                 if(anchorBitPairs < 32) {
1576                         clearMask >>= ((32-anchorBitPairs) << 1);
1577                         useMask = true;
1578                 }
1579                 int nsInAnchor = 0;
1580                 uint32_t nPoss = 0;
1581                 int nPos1 = -1;
1582                 int nPos2 = -1;
1583                 int nPos3 = -1;
1584                 uint32_t skipLeftToRights = 0;
1585                 uint32_t skipRightToLefts = 0;
1586                 // Construct the 'anchor' 64-bit buffer so that it holds all of
1587                 // the first 'anchorBitPairs' bit pairs of the query.
1588                 for(uint32_t i = 0; i < anchorBitPairs; i++) {
1589                         int c = (int)qry[i]; // next query character
1590                         int r = (int)ref[halfway - begin + i]; // next reference character
1591                         if(r & 4) {
1592                                 r = 0;
1593                                 // The right-to-left direction absorbs the candidate
1594                                 // alignment based at 'halfway'; so skipLeftToRights is
1595                                 // i, not i+1
1596                                 skipLeftToRights = max(skipLeftToRights, i);
1597                                 skipRightToLefts = max(skipRightToLefts, anchorBitPairs - i);
1598                         }
1599                         assert_leq(c, 4);
1600                         assert_lt(r, 4);
1601                         // Special case: query has an 'N'
1602                         if(c == 4) {
1603                                 if(++nsInAnchor > 3) {
1604                                         // More than two 'N's in the anchor region; can't
1605                                         // possibly have a 2-mismatch hit anywhere
1606                                         assert_eq(r2.size(), ranges.size() - rangesInitSz);
1607                                         return;   // can't match if query has Ns
1608                                 } else if(nsInAnchor == 3) {
1609                                         nPos3 = (int)i;
1610                                         nPoss++;
1611                                 } else if(nsInAnchor == 2) {
1612                                         nPos2 = (int)i;
1613                                         nPoss++;
1614                                 } else {
1615                                         assert_eq(1, nsInAnchor);
1616                                         nPos1 = (int)i;
1617                                         nPoss++;
1618                                 }
1619                                 // Make it look like an 'A' in the anchor
1620                                 c = 0;
1621                                 diffMask = (diffMask << 2llu) | 1llu;
1622                         } else {
1623                                 diffMask <<= 2llu;
1624                         }
1625                         anchor  = ((anchor  << 2llu) | c);
1626                         buffw = ((buffw << 2llu) | r);
1627                 }
1628                 assert_leq(nPoss, 3);
1629                 // Check whether read is disqualified by Ns outside of the anchor
1630                 // region
1631                 for(uint32_t i = anchorBitPairs; i < qlen; i++) {
1632                         if((int)qry[i] == 4) {
1633                                 if(++nsInAnchor > 3) {
1634                                         assert_eq(r2.size(), ranges.size() - rangesInitSz);
1635                                         return; // can't match if query has Ns
1636                                 }
1637                         }
1638                 }
1639                 uint64_t bufbw = buffw;
1640                 // We're moving the right-hand edge of the anchor along until
1641                 // it's 'anchorOverhang' chars from the end of the target region.
1642                 // Note that we're not making a 3'/5' distinction here; if we
1643                 // were, we might need to make the 'anchorOverhang' adjustment on
1644                 // the left end of the range rather than the right.
1645                 bool hi = false;
1646                 uint32_t riHi  = halfway;
1647                 uint32_t rirHi = halfway - begin;
1648                 uint32_t rirHiAnchor = rirHi + anchorBitPairs - 1;
1649                 uint32_t riLo  = halfway + 1;
1650                 uint32_t rirLo = halfway - begin + 1;
1651                 for(uint32_t i = 1; i <= lim + 1; i++) {
1652                         int r;       // new reference char
1653                         uint64_t diff;
1654                         assert_lt(skipLeftToRights, qlen);
1655                         assert_leq(skipRightToLefts, qlen);
1656                         if(hi) {
1657                                 hi = false;
1658                                 // Moving left-to-right
1659                                 riHi++;
1660                                 rirHi++;
1661                                 rirHiAnchor++;
1662                                 r = (int)ref[rirHiAnchor];
1663                                 if(r & 4) {
1664                                         r = 0;
1665                                         skipLeftToRights = anchorBitPairs;
1666                                 }
1667                                 assert_lt(r, 4);
1668                                 // Bring in new base pair at the least significant
1669                                 // position
1670                                 buffw = ((buffw << 2llu) | r);
1671                                 if(useMask) buffw &= clearMask;
1672                                 if(skipLeftToRights > 0) {
1673                                         skipLeftToRights--;
1674                                         continue;
1675                                 }
1676                                 diff = (buffw ^ anchor) | diffMask;
1677                         } else {
1678                                 hi = true;
1679                                 // Moving right-to-left
1680                                 riLo--;
1681                                 rirLo--;
1682                                 r = (int)ref[rirLo];
1683                                 if(r & 4) {
1684                                         r = 0;
1685                                         skipRightToLefts = qlen;
1686                                 }
1687                                 assert_lt(r, 4);
1688                                 if(i >= 2) {
1689                                         bufbw >>= 2llu;
1690                                         // Bring in new base pair at the most significant
1691                                         // position
1692                                         bufbw |= ((uint64_t)r << lhsShift);
1693                                 }
1694                                 if(skipRightToLefts > 0) {
1695                                         skipRightToLefts--;
1696                                         continue;
1697                                 }
1698                                 diff = (bufbw ^ anchor) | diffMask;
1699                         }
1700                         if((diff & 0xffff000000000000llu) &&
1701                            (diff & 0x0000ffff00000000llu) &&
1702                            (diff & 0x00000000ffff0000llu) &&
1703                            (diff & 0x000000000000ffffllu)) continue;
1704                         uint32_t ri  = hi ? riLo  : riHi;
1705                         uint32_t rir = hi ? rirLo : rirHi;
1706                         // Could use pop count
1707                         uint8_t *diff8 = reinterpret_cast<uint8_t*>(&diff);
1708                         // As a first cut, see if there are too many mismatches in
1709                         // the first and last parts of the anchor
1710                         uint32_t diffs = u8toMms[(int)diff8[0]] + u8toMms[(int)diff8[7]];
1711                         if(diffs > 3) continue;
1712                         diffs += u8toMms[(int)diff8[1]] +
1713                                          u8toMms[(int)diff8[2]] +
1714                                          u8toMms[(int)diff8[3]] +
1715                                          u8toMms[(int)diff8[4]] +
1716                                          u8toMms[(int)diff8[5]] +
1717                                          u8toMms[(int)diff8[6]];
1718                         uint32_t mmpos1 = 0xffffffff;
1719                         int refc1 = -1;
1720                         uint32_t mmpos2 = 0xffffffff;
1721                         int refc2 = -1;
1722                         uint32_t mmpos3 = 0xffffffff;
1723                         int refc3 = -1;
1724                         if(diffs > 3) {
1725                                 // Too many differences
1726                                 continue;
1727                         } else if(nPoss > 1 && diffs == nPoss) {
1728                                 // There was one difference, but there was also one N,
1729                                 // so we already know where the difference is
1730                                 mmpos1 = nPos1;
1731                                 refc1 = "ACGT"[(int)ref[rir + nPos1]];
1732                                 if(nPoss > 1) {
1733                                         mmpos2 = nPos2;
1734                                         refc2 = "ACGT"[(int)ref[rir + nPos2]];
1735                                         if(nPoss > 2) {
1736                                                 mmpos3 = nPos3;
1737                                                 refc3 = "ACGT"[(int)ref[rir + nPos3]];
1738                                         }
1739                                 }
1740                         } else if(diffs > 0) {
1741                                 // Figure out which position mismatched
1742                                 uint64_t diff2 = diff;
1743                                 mmpos1 = 31;
1744                                 if((diff & 0xffffffffllu) == 0) { diff >>= 32llu; mmpos1 -= 16; }
1745                                 assert_neq(0, diff);
1746                                 if((diff & 0xffffllu) == 0)     { diff >>= 16llu; mmpos1 -=  8; }
1747                                 assert_neq(0, diff);
1748                                 if((diff & 0xffllu) == 0)       { diff >>= 8llu;  mmpos1 -=  4; }
1749                                 assert_neq(0, diff);
1750                                 if((diff & 0xfllu) == 0)        { diff >>= 4llu;  mmpos1 -=  2; }
1751                                 assert_neq(0, diff);
1752                                 if((diff & 0x3llu) == 0)        { mmpos1--; }
1753                                 assert_neq(0, diff);
1754                                 assert_geq(mmpos1, 0);
1755                                 assert_lt(mmpos1, 32);
1756                                 mmpos1 -= anchorCushion;
1757                                 refc1 = "ACGT"[(int)ref[rir + mmpos1]];
1758                                 if(diffs > 1) {
1759                                         // Figure out the second mismatched position
1760                                         diff2 &= ~(0xc000000000000000llu >> (uint64_t)((mmpos1 + anchorCushion) << 1));
1761                                         uint64_t diff3 = diff2;
1762                                         mmpos2 = 31;
1763                                         if((diff2 & 0xffffffffllu) == 0) { diff2 >>= 32llu; mmpos2 -= 16; }
1764                                         assert_neq(0, diff2);
1765                                         if((diff2 & 0xffffllu) == 0)     { diff2 >>= 16llu; mmpos2 -=  8; }
1766                                         assert_neq(0, diff2);
1767                                         if((diff2 & 0xffllu) == 0)       { diff2 >>= 8llu;  mmpos2 -=  4; }
1768                                         assert_neq(0, diff2);
1769                                         if((diff2 & 0xfllu) == 0)        { diff2 >>= 4llu;  mmpos2 -=  2; }
1770                                         assert_neq(0, diff2);
1771                                         if((diff2 & 0x3llu) == 0)        { mmpos2--; }
1772                                         assert_neq(0, diff2);
1773                                         mmpos2 -= anchorCushion;
1774                                         assert_geq(mmpos2, 0);
1775                                         assert_lt(mmpos2, 32);
1776                                         assert_neq(mmpos1, mmpos2);
1777                                         refc2 = "ACGT"[(int)ref[rir + mmpos2]];
1778                                         uint32_t mmpos2orig = mmpos2;
1779                                         if(mmpos2 < mmpos1) {
1780                                                 uint32_t mmtmp = mmpos1;
1781                                                 mmpos1 = mmpos2;
1782                                                 mmpos2 = mmtmp;
1783                                                 int refctmp = refc1;
1784                                                 refc1 = refc2;
1785                                                 refc2 = refctmp;
1786                                         }
1787                                         assert_lt(mmpos1, mmpos2);
1788                                         if(diffs > 2) {
1789                                                 // Figure out the second mismatched position
1790                                                 diff3 &= ~(0xc000000000000000llu >> (uint64_t)((mmpos2orig + anchorCushion) << 1));
1791                                                 mmpos3 = 31;
1792                                                 if((diff3 & 0xffffffffllu) == 0) { diff3 >>= 32llu; mmpos3 -= 16; }
1793                                                 assert_neq(0, diff3);
1794                                                 if((diff3 & 0xffffllu) == 0)     { diff3 >>= 16llu; mmpos3 -=  8; }
1795                                                 assert_neq(0, diff3);
1796                                                 if((diff3 & 0xffllu) == 0)       { diff3 >>= 8llu;  mmpos3 -=  4; }
1797                                                 assert_neq(0, diff3);
1798                                                 if((diff3 & 0xfllu) == 0)        { diff3 >>= 4llu;  mmpos3 -=  2; }
1799                                                 assert_neq(0, diff3);
1800                                                 if((diff3 & 0x3llu) == 0)        { mmpos3--; }
1801                                                 assert_neq(0, diff3);
1802                                                 mmpos3 -= anchorCushion;
1803                                                 assert_geq(mmpos3, 0);
1804                                                 assert_lt(mmpos3, 32);
1805                                                 assert_neq(mmpos1, mmpos3);
1806                                                 assert_neq(mmpos2, mmpos3);
1807                                                 refc3 = "ACGT"[(int)ref[rir + mmpos3]];
1808                                                 if(mmpos3 < mmpos1) {
1809                                                         uint32_t mmtmp = mmpos1;
1810                                                         mmpos1 = mmpos3;
1811                                                         mmpos3 = mmpos2;
1812                                                         mmpos2 = mmtmp;
1813                                                         int refctmp = refc1;
1814                                                         refc1 = refc3;
1815                                                         refc3 = refc2;
1816                                                         refc2 = refctmp;
1817                                                 } else if(mmpos3 < mmpos2) {
1818                                                         uint32_t mmtmp = mmpos2;
1819                                                         mmpos2 = mmpos3;
1820                                                         mmpos3 = mmtmp;
1821                                                         int refctmp = refc2;
1822                                                         refc2 = refc3;
1823                                                         refc3 = refctmp;
1824                                                 }
1825                                                 assert_lt(mmpos1, mmpos2);
1826                                                 assert_lt(mmpos2, mmpos3);
1827                                         }
1828                                 }
1829                         }
1830                         // Now extend the anchor into a longer alignment
1831                         bool foundHit = true;
1832                         if(anchorOverhang > 0) {
1833                                 assert_leq(ri + anchorBitPairs + anchorOverhang, end);
1834                                 bool skipCandidate = false;
1835                                 for(uint32_t j = 0; j < anchorOverhang; j++) {
1836                                         int rc = (int)ref[rir + 32 + j];
1837                                         if(rc == 4) {
1838                                                 // Oops, encountered an N in the reference in
1839                                                 // the overhang portion of the candidate
1840                                                 // alignment
1841                                                 // (Note that we inverted hi earlier)
1842                                                 if(hi) {
1843                                                         // Right-to-left
1844                                                         skipRightToLefts = anchorOverhang - j - 1;
1845                                                 } else {
1846                                                         // Left-to-right
1847                                                         skipLeftToRights = anchorBitPairs + j;
1848                                                 }
1849                                                 skipCandidate = true; // Skip this candidate
1850                                                 break;
1851                                         }
1852                                         if((int)qry[32 + j] != rc) {
1853                                                 if(++diffs > 3) {
1854                                                         foundHit = false;
1855                                                         break;
1856                                                 } else if(diffs == 3) {
1857                                                         mmpos3 = 32 + j;
1858                                                         refc3 = "ACGT"[(int)ref[rir + 32 + j]];
1859                                                 } else if(diffs == 2) {
1860                                                         mmpos2 = 32 + j;
1861                                                         refc2 = "ACGT"[(int)ref[rir + 32 + j]];
1862                                                 } else {
1863                                                         assert_eq(1, diffs);
1864                                                         mmpos1 = 32 + j;
1865                                                         refc1 = "ACGT"[(int)ref[rir + 32 + j]];
1866                                                 }
1867                                         }
1868                                 }
1869                                 if(skipCandidate) continue;
1870                         }
1871                         if(!foundHit) continue;
1872                         if(pairs != NULL) {
1873                                 TU64Pair p;
1874                                 if(ri < aoff) {
1875                                         // By convention, the upstream mate's
1876                                         // coordinates go in the 'first' field
1877                                         p.first  = ((uint64_t)tidx << 32) | (uint64_t)ri;
1878                                         p.second = ((uint64_t)tidx << 32) | (uint64_t)aoff;
1879                                 } else {
1880                                         p.second = ((uint64_t)tidx << 32) | (uint64_t)ri;
1881                                         p.first  = ((uint64_t)tidx << 32) | (uint64_t)aoff;
1882                                 }
1883                                 if(pairs->find(p) != pairs->end()) {
1884                                         // We already found this hit!  Continue.
1885                                         ASSERT_ONLY(duplicates++);
1886                                         ASSERT_ONLY(r2i++);
1887                                         continue;
1888                                 } else {
1889                                         // Record this hit
1890                                         pairs->insert(p);
1891                                 }
1892                         }
1893                         assert_leq(diffs, 3);
1894                         assert_lt(r2i, r2.size());
1895                         assert_eq(re2[r2i], ri);
1896                         ranges.resize(ranges.size()+1);
1897                         Range& range = ranges.back();
1898                         assert_eq(diffs, r2[r2i].numMms);
1899                         range.stratum = diffs;
1900                         range.numMms = diffs;
1901                         assert_eq(0, range.mms.size());
1902                         assert_eq(0, range.refcs.size());
1903                         if(diffs > 0) {
1904                                 assert_neq(mmpos1, 0xffffffff);
1905                                 assert_eq(mmpos1, r2[r2i].mms[0]);
1906                                 assert_neq(-1, refc1);
1907                                 assert_eq(refc1, r2[r2i].refcs[0]);
1908                                 range.mms.push_back(mmpos1);
1909                                 range.refcs.push_back(refc1);
1910                                 if(diffs > 1) {
1911                                         assert_neq(mmpos2, 0xffffffff);
1912                                         assert_eq(mmpos2, r2[r2i].mms[1]);
1913                                         assert_neq(-1, refc2);
1914                                         assert_eq(refc2, r2[r2i].refcs[1]);
1915                                         range.mms.push_back(mmpos2);
1916                                         range.refcs.push_back(refc2);
1917                                         if(diffs > 2) {
1918                                                 assert_neq(mmpos3, 0xffffffff);
1919                                                 assert_eq(mmpos3, r2[r2i].mms[2]);
1920                                                 assert_neq(-1, refc3);
1921                                                 assert_eq(refc3, r2[r2i].refcs[2]);
1922                                                 range.mms.push_back(mmpos3);
1923                                                 range.refcs.push_back(refc3);
1924                                         }
1925                                 }
1926                         }
1927                         ASSERT_ONLY(r2i++);
1928                         results.push_back(ri);
1929                         if(--numToFind == 0) return;
1930                 }
1931                 assert_leq(duplicates, r2.size());
1932                 assert_geq(r2.size() - duplicates, ranges.size() - rangesInitSz);
1933                 return; // no hit
1934         }
1935 };
1936
1937 /**
1938  * Concrete RefAligner for finding nearby 1-mismatch hits given an
1939  * anchor hit.
1940  */
1941 template<typename TStr>
1942 class Seed0RefAligner : public RefAligner<TStr> {
1943
1944         typedef seqan::String<seqan::Dna5> TDna5Str;
1945         typedef seqan::String<char> TCharStr;
1946         typedef std::vector<Range> TRangeVec;
1947         typedef std::vector<uint32_t> TU32Vec;
1948         typedef std::pair<uint64_t, uint64_t> TU64Pair;
1949         typedef std::set<TU64Pair> TSetPairs;
1950
1951 public:
1952
1953         Seed0RefAligner(bool color, bool verbose, bool quiet, uint32_t seedLen, uint32_t qualMax, bool maqPenalty) :
1954                 RefAligner<TStr>(color, verbose, quiet, seedLen, qualMax, maqPenalty) { }
1955
1956         virtual ~Seed0RefAligner() { }
1957
1958 protected:
1959         /**
1960          * This schematic shows the roles played by the begin, qbegin, end,
1961          * qend, halfway, slen, qlen, and lim variables:
1962          *
1963          * seedOnLeft == true:
1964          *
1965          * |<                   lim                   >|<     qlen       >|
1966          *  --------------------------------------------------------------
1967          * |                     | slen | qlen-slen |  | slen | qlen-slen |
1968          *  --------------------------------------------------------------
1969          * ^                     ^                     ^                  ^
1970          * begin & qbegin     halfway                qend               end
1971          *
1972          * seedOnLeft == false:
1973          *
1974          * |<     qlen       >|<                   lim                   >|
1975          *  --------------------------------------------------------------
1976          * | qlen-slen | slen |  | qlen-slen | slen |                     |
1977          *  --------------------------------------------------------------
1978          * ^                  ^                     ^                     ^
1979          * begin            qbegin             halfway           qend & end
1980          *
1981          * Note that, for seeds longer than 32 base-pairs, the seed is
1982          * further subdivided into a 32-bit anchor and a seed overhang of
1983          * length > 0.
1984          */
1985         void naiveFind(uint32_t numToFind,
1986                                    uint32_t tidx,
1987                                    uint8_t* ref,
1988                                    const TDna5Str& qry,
1989                                    const TCharStr& quals,
1990                                    uint32_t begin,
1991                                    uint32_t end,
1992                                    TRangeVec& ranges,
1993                                    TU32Vec& results,
1994                                    TSetPairs* pairs,
1995                                    uint32_t aoff,
1996                                    bool seedOnLeft) const
1997         {
1998                 assert_gt(numToFind, 0);
1999                 assert_gt(end, begin);
2000                 const uint32_t qlen = seqan::length(qry);
2001                 assert_gt(qlen, 0);
2002                 assert_geq(end - begin, qlen); // caller should have checked this
2003                 assert_gt(this->seedLen_, 0);
2004                 const uint32_t slen = min(qlen, this->seedLen_);
2005                 uint32_t qend = end;
2006                 uint32_t qbegin = begin;
2007                 // If the seed is on the left-hand side of the alignment, then
2008                 // leave a gap at the right-hand side of the interval;
2009                 // otherwise, do the opposite
2010                 if(seedOnLeft) {
2011                         // Leave gap on right-hand side of the interval
2012                         qend -= qlen;
2013                 } else {
2014                         // Leave gap on left-hand side of the interval
2015                         qbegin += qlen;
2016                 }
2017                 // lim = number of alignments to try
2018                 const uint32_t lim = qend - qbegin;
2019                 // halfway = position in the reference to start at (and then
2020                 // we work our way out to the right and to the left).
2021                 const uint32_t halfway = qbegin + (lim >> 1);
2022                 // Vectors for holding edit information
2023                 std::vector<uint32_t> nonSeedMms;
2024                 std::vector<uint8_t> nonSeedRefcs;
2025                 bool hi = false;
2026                 for(uint32_t i = 1; i <= lim+1; i++) {
2027                         uint32_t ri;  // leftmost position in candidate alignment
2028                         uint32_t rir; // same, minus begin; for indexing into ref[]
2029                         if(hi) {
2030                                 ri = halfway + (i >> 1); rir = ri - begin;
2031                                 assert_leq(ri, qend);
2032                         } else {
2033                                 ri = halfway - (i >> 1); rir = ri - begin;
2034                                 assert_geq(ri, begin);
2035                         }
2036                         hi = !hi;
2037                         // Do the naive comparison
2038                         bool match = true;
2039                         int mms = 0;
2040                         unsigned int ham = 0;
2041                         nonSeedMms.clear();
2042                         nonSeedRefcs.clear();
2043                         // Walk through each position of the alignment
2044                         for(uint32_t jj = 0; jj < qlen; jj++) {
2045                                 uint32_t j = jj;
2046                                 if(!seedOnLeft) {
2047                                         // If seed is on the right, scan right-to-left
2048                                         j = qlen - jj - 1;
2049                                 } else {
2050                                         // Go left-to-right
2051                                 }
2052                                 uint32_t rirj = rir + j;
2053                                 if(!seedOnLeft) {
2054                                         assert_geq(rir, jj);
2055                                         rirj = rir - jj - 1;
2056                                 }
2057 #if 0
2058                                 // Count Ns in the reference as mismatches
2059                                 const int q = (int)qry[j];
2060                                 const int r = (int)ref[rirj];
2061                                 assert_leq(q, 4);
2062                                 assert_leq(r, 4);
2063                                 if(q == 4 || r == 4 || q != r) {
2064 #else
2065                                 // Disallow alignments that involve an N in the
2066                                 // reference
2067                                 const int r = (int)ref[rirj];
2068                                 if(r & 4) {
2069                                         // N in reference; bail on this alignment
2070                                         match = false;
2071                                         break;
2072                                 }
2073                                 const int q = (int)qry[j];
2074                                 assert_leq(q, 4);
2075                                 assert_lt(r, 4);
2076                                 if(q != r) {
2077 #endif
2078                                         // Mismatch!
2079                                         if(jj < slen) {
2080                                                 // More than one mismatch in the anchor; reject
2081                                                 match = false;
2082                                                 break;
2083                                         }
2084                                         uint8_t qual = phredCharToPhredQual(quals[j]);
2085                                         ham += mmPenalty(this->maqPenalty_, qual);
2086                                         if(ham > this->qualMax_) {
2087                                                 // Exceeded quality ceiling; reject
2088                                                 match = false;
2089                                                 break;
2090                                         } else {
2091                                                 // Legal mismatch outside of the anchor; record it
2092                                                 mms++;
2093                                                 nonSeedMms.push_back(j);
2094                                                 assert_leq(nonSeedMms.size(), (size_t)mms);
2095                                                 nonSeedRefcs.push_back("ACGTN"[r]);
2096                                         }
2097                                 }
2098                         }
2099                         if(match) {
2100                                 ranges.resize(ranges.size()+1);
2101                                 Range& range = ranges.back();
2102                                 range.stratum = 0;
2103                                 range.numMms = mms;
2104                                 assert_eq(0, range.mms.size());
2105                                 assert_eq(0, range.refcs.size());
2106                                 if(mms >= 1) {
2107                                         // Be careful to add edits in left-to-right order
2108                                         // with respect to the read/alignment
2109                                         const size_t nonSeedMmsSz = nonSeedMms.size();
2110                                         if(nonSeedMmsSz > 0) {
2111                                                 if(seedOnLeft) {
2112                                                         for(size_t k = 0; k < nonSeedMmsSz; k++) {
2113                                                                 range.mms.push_back(nonSeedMms[k]);
2114                                                                 range.refcs.push_back(nonSeedRefcs[k]);
2115                                                         }
2116                                                 } else {
2117                                                         for(size_t k = nonSeedMmsSz; k > 0; k--) {
2118                                                                 range.mms.push_back(nonSeedMms[k-1]);
2119                                                                 range.refcs.push_back(nonSeedRefcs[k-1]);
2120                                                         }
2121                                                 }
2122                                         }
2123                                 }
2124                                 assert_eq((size_t)mms, range.mms.size());
2125                                 assert_eq((size_t)mms, range.refcs.size());
2126                                 if(seedOnLeft) {
2127                                         results.push_back(ri);
2128                                 } else {
2129                                         results.push_back(ri - qlen);
2130                                 }
2131                         }
2132                 }
2133                 return;
2134         }
2135
2136         /**
2137          * This schematic shows the roles played by the begin, qbegin, end,
2138          * qend, halfway, slen, qlen, and lim variables:
2139          *
2140          * seedOnLeft == true:
2141          *
2142          * |<                   lim                   >|<     qlen       >|
2143          *  --------------------------------------------------------------
2144          * |                     | slen | qlen-slen |  | slen | qlen-slen |
2145          *  --------------------------------------------------------------
2146          * ^                     ^                     ^                  ^
2147          * begin & qbegin     halfway                qend               end
2148          *
2149          * seedOnLeft == false:
2150          *
2151          *             |<                   lim                   >|
2152          *  --------------------------------------------------------------
2153          * | qlen-slen |         | qlen-slen | slen |              | slen |
2154          *  --------------------------------------------------------------
2155          * ^           ^                     ^                     ^      ^
2156          * begin       qbegin             halfway                qend   end
2157          *
2158          * Note that, for seeds longer than 32 base-pairs, the seed is
2159          * further subdivided into a 32-bit anchor and a seed overhang of
2160          * length > 0.
2161          */
2162         virtual void anchor64Find(uint32_t numToFind,
2163                                         uint32_t tidx,
2164                                         uint8_t* ref,
2165                                         const TDna5Str& qry,
2166                                         const TCharStr& quals,
2167                                         uint32_t begin,
2168                                         uint32_t end,
2169                                         TRangeVec& ranges,
2170                                         TU32Vec& results,
2171                                         TSetPairs* pairs = NULL,
2172                                         uint32_t aoff = 0xffffffff,
2173                                         bool seedOnLeft = false) const
2174         {
2175                 assert_gt(numToFind, 0);
2176                 ASSERT_ONLY(const uint32_t rangesInitSz = ranges.size());
2177                 ASSERT_ONLY(uint32_t duplicates = 0);
2178                 ASSERT_ONLY(uint32_t r2i = 0);
2179                 const uint32_t qlen = seqan::length(qry);
2180                 assert_gt(qlen, 0);
2181                 assert_gt(end, begin);
2182                 assert_geq(end - begin, qlen); // caller should have checked this
2183                 assert_gt(this->seedLen_, 0);
2184                 uint32_t slen = min(qlen, this->seedLen_);
2185 #ifndef NDEBUG
2186                 // Get results from the naive matcher for sanity-checking
2187                 TRangeVec r2; TU32Vec re2;
2188                 naiveFind(numToFind, tidx, ref, qry, quals, begin, end, r2,
2189                                   re2, pairs, aoff, seedOnLeft);
2190 #endif
2191                 const uint32_t anchorBitPairs = min<int>(slen, 32);
2192                 const int lhsShift = ((anchorBitPairs - 1) << 1);
2193                 ASSERT_ONLY(const uint32_t anchorCushion  = 32 - anchorBitPairs);
2194                 // seedAnchorOverhang = # seed bases not included in the anchor
2195                 const uint32_t seedAnchorOverhang = (slen <= anchorBitPairs ? 0 : (slen - anchorBitPairs));
2196                 // seedAnchorOverhang = # seed bases not included in the anchor
2197                 const uint32_t readSeedOverhang = (slen == qlen ? 0 : (qlen - slen));
2198                 assert(anchorCushion == 0 || seedAnchorOverhang == 0);
2199                 assert_eq(qlen, readSeedOverhang + slen);
2200                 uint32_t qend = end;
2201                 uint32_t qbegin = begin;
2202                 if(seedOnLeft) {
2203                         // Leave read-sized gap on right-hand side of the interval
2204                         qend -= qlen;
2205                 } else {
2206                         // Leave seed-sized gap on right-hand side and
2207                         // non-seed-sized gap on the left-hand side
2208                         qbegin += readSeedOverhang;
2209                         qend -= slen;
2210                 }
2211                 // lim = # possible alignments in the range
2212                 const uint32_t lim = qend - qbegin;
2213                 // halfway = point on the genome to radiate out from
2214                 const uint32_t halfway = qbegin + (lim >> 1);
2215                 uint64_t anchor = 0llu;
2216                 uint64_t buffw = 0llu;  // rotating ref sequence buffer
2217                 // Set up a mask that we'll apply to the two bufs every round
2218                 // to discard bits that were rotated out of the anchor area
2219                 uint64_t clearMask = 0xffffffffffffffffllu;
2220                 bool useMask = false;
2221                 if(anchorBitPairs < 32) {
2222                         clearMask >>= ((32-anchorBitPairs) << 1);
2223                         useMask = true;
2224                 }
2225                 uint32_t skipLeftToRights = 0;
2226                 uint32_t skipRightToLefts = 0;
2227                 const uint32_t halfwayRi = halfway - begin;
2228                 // Construct the 'anchor' 64-bit buffer so that it holds all of
2229                 // the first 'anchorBitPairs' bit pairs of the query.
2230                 for(uint32_t ii = 0; ii < anchorBitPairs; ii++) {
2231                         uint32_t i = ii;
2232                         if(!seedOnLeft) {
2233                                 // Fill in the anchor using characters from the right-
2234                                 // hand side of the query (but take the characters in
2235                                 // left-to-right order)
2236                                 i = qlen - slen + ii;
2237                         }
2238                         int c = (int)qry[i]; // next query character
2239                         int r = (int)ref[halfwayRi + ii]; // next reference character
2240                         if(r & 4) {
2241                                 // The reference character is an N; to mimic the
2242                                 // behavior of BW alignment, we have to skip all
2243                                 // alignments that involve an N in the reference.  Set
2244                                 // the skip* variables accordingly.
2245                                 r = 0;
2246                                 uint32_t lrSkips = ii;
2247                                 uint32_t rlSkips = qlen - ii;
2248                                 if(!seedOnLeft && readSeedOverhang) {
2249                                         lrSkips += readSeedOverhang;
2250                                         assert_geq(rlSkips, readSeedOverhang);
2251                                         rlSkips -= readSeedOverhang;
2252                                 }
2253                                 // The right-to-left direction absorbs the candidate
2254                                 // alignment based at 'halfway'
2255                                 skipLeftToRights = max(skipLeftToRights, lrSkips);
2256                                 skipRightToLefts = max(skipRightToLefts, rlSkips);
2257                                 assert_leq(skipLeftToRights, qlen);
2258                                 assert_leq(skipRightToLefts, qlen);
2259                         }
2260                         assert_leq(c, 4);
2261                         assert_lt(r, 4);
2262                         // Special case: query has an 'N'
2263                         if(c == 4) {
2264                                 // One or more 'N's in the anchor region; can't
2265                                 // possibly have a 0-mismatch hit anywhere
2266                                 assert_eq(r2.size(), ranges.size() - rangesInitSz);
2267                                 return;   // can't match if query has Ns
2268                         }
2269                         anchor = ((anchor << 2llu) | c);
2270                         buffw = ((buffw << 2llu) | r);
2271                 }
2272                 // Check whether read is disqualified by Ns inside the seed
2273                 // region but outside the anchor region
2274                 if(seedAnchorOverhang) {
2275                         assert_lt(anchorBitPairs, slen);
2276                         for(uint32_t ii = anchorBitPairs; ii < slen; ii++) {
2277                                 uint32_t i = ii;
2278                                 if(!seedOnLeft) {
2279                                         i = qlen - slen + ii;
2280                                 }
2281                                 if((int)qry[i] == 4) {
2282                                         assert_eq(r2.size(), ranges.size() - rangesInitSz);
2283                                         return; // can't match if query has Ns
2284                                 }
2285                         }
2286                 } else {
2287                         assert_eq(anchorBitPairs, slen);
2288                 }
2289                 uint64_t bufbw = buffw;
2290                 // Slide the anchor out in either direction, alternating
2291                 // between right-to-left and left-to-right shifts, until all of
2292                 // the positions from qbegin to qend have been covered.
2293                 bool hi = false;
2294                 uint32_t riHi  = halfway;
2295                 uint32_t rirHi = halfway - begin;
2296                 uint32_t rirHiAnchor = rirHi + anchorBitPairs - 1;
2297                 uint32_t riLo  = halfway + 1;
2298                 uint32_t rirLo = halfway - begin + 1;
2299                 uint32_t lrSkips = anchorBitPairs;
2300                 uint32_t rlSkips = qlen;
2301                 if(!seedOnLeft && readSeedOverhang) {
2302                         lrSkips += readSeedOverhang;
2303                         assert_geq(rlSkips, readSeedOverhang);
2304                         rlSkips -= readSeedOverhang;
2305                 }
2306                 for(uint32_t i = 1; i <= lim + 1; i++) {
2307                         int r;       // new reference char
2308                         uint64_t diff;
2309                         assert_leq(skipLeftToRights, qlen);
2310                         assert_leq(skipRightToLefts, qlen);
2311                         if(hi) {
2312                                 hi = false;
2313                                 // Moving left-to-right
2314                                 riHi++;
2315                                 rirHi++;
2316                                 rirHiAnchor++;
2317                                 r = (int)ref[rirHiAnchor];
2318                                 if(r & 4) {
2319                                         r = 0;
2320                                         skipLeftToRights = lrSkips;
2321                                 }
2322                                 assert_lt(r, 4);
2323                                 // Bring in new base pair at the least significant
2324                                 // position
2325                                 buffw = ((buffw << 2llu) | r);
2326                                 if(useMask) buffw &= clearMask;
2327                                 if(skipLeftToRights > 0) {
2328                                         skipLeftToRights--;
2329                                         continue;
2330                                 }
2331                                 diff = buffw ^ anchor;
2332                         } else {
2333                                 hi = true;
2334                                 // Moving right-to-left
2335                                 riLo--;
2336                                 rirLo--;
2337                                 r = (int)ref[rirLo];
2338                                 if(r & 4) {
2339                                         r = 0;
2340                                         skipRightToLefts = rlSkips;
2341                                 }
2342                                 assert_lt(r, 4);
2343                                 if(i >= 2) {
2344                                         bufbw >>= 2llu;
2345                                         // Bring in new base pair at the most significant
2346                                         // position
2347                                         bufbw |= ((uint64_t)r << lhsShift);
2348                                 }
2349                                 if(skipRightToLefts > 0) {
2350                                         skipRightToLefts--;
2351                                         continue;
2352                                 }
2353                                 diff = bufbw ^ anchor;
2354                         }
2355                         if(diff) continue;
2356                         uint32_t ri  = hi ? riLo  : riHi;
2357                         uint32_t rir = hi ? rirLo : rirHi;
2358                         unsigned int ham = 0;
2359                         // If the seed is longer than the anchor, then scan the
2360                         // rest of the seed characters
2361                         bool foundHit = true;
2362                         if(seedAnchorOverhang) {
2363                                 for(uint32_t j = 0; j < seedAnchorOverhang; j++) {
2364                                         int rc = (int)ref[rir + anchorBitPairs + j];
2365                                         if(rc == 4) {
2366                                                 // Oops, encountered an N in the reference in
2367                                                 // the overhang portion of the candidate
2368                                                 // alignment
2369                                                 // (Note that we inverted hi earlier)
2370                                                 if(hi) {
2371                                                         // Right-to-left
2372                                                         // Skip out of the seedAnchorOverhang
2373                                                         assert_eq(0, skipRightToLefts);
2374                                                         skipRightToLefts = seedAnchorOverhang - j - 1;
2375                                                         if(seedOnLeft) {
2376                                                                 // ...and skip out of the rest of the read
2377                                                                 skipRightToLefts += readSeedOverhang;
2378                                                         }
2379                                                 } else {
2380                                                         // Left-to-right
2381                                                         // Skip out of the seedAnchorOverhang
2382                                                         assert_eq(0, skipLeftToRights);
2383                                                         skipLeftToRights = anchorBitPairs + j;
2384                                                         if(!seedOnLeft) {
2385                                                                 // ...and skip out of the rest of the read
2386                                                                 skipLeftToRights += readSeedOverhang;
2387                                                         }
2388                                                 }
2389                                                 foundHit = false; // Skip this candidate
2390                                                 break;
2391                                         }
2392                                         uint32_t qoff = anchorBitPairs + j;
2393                                         if(!seedOnLeft) {
2394                                                 qoff += readSeedOverhang;
2395                                         }
2396                                         assert_lt(qoff, qlen);
2397                                         if((int)qry[qoff] != rc) {
2398                                                 foundHit = false;
2399                                                 break;
2400                                         }
2401                                 }
2402                                 if(!foundHit) continue;
2403                         }
2404                         // If the read is longer than the seed, then scan the rest
2405                         // of the read characters; mismatches no longer count
2406                         // toward the stratum or the 1-mm limit.
2407                         // Vectors for holding edit information
2408                         std::vector<uint32_t> nonSeedMms;
2409                         std::vector<uint8_t> nonSeedRefcs;
2410                         int mms = 0; // start counting total mismatches
2411                         if((qlen - slen) > 0) {
2412                                 // Going left-to-right
2413                                 for(uint32_t j = 0; j < readSeedOverhang; j++) {
2414                                         uint32_t roff = rir + slen + j;
2415                                         uint32_t qoff = slen + j;
2416                                         if(!seedOnLeft) {
2417                                                 assert_geq(roff, qlen);
2418                                                 roff -= qlen;
2419                                                 qoff = j;
2420                                         }
2421                                         int rc = (int)ref[roff];
2422                                         if(rc == 4) {
2423                                                 // Oops, encountered an N in the reference in
2424                                                 // the overhang portion of the candidate
2425                                                 // alignment
2426                                                 // (Note that we inverted hi earlier)
2427                                                 if(hi) {
2428                                                         // Right-to-left
2429                                                         // Skip what's left of the readSeedOverhang
2430                                                         skipRightToLefts = readSeedOverhang - j - 1;
2431                                                         if(!seedOnLeft) {
2432                                                                 // ...and skip the seed if it's on the right
2433                                                                 skipRightToLefts += slen;
2434                                                         }
2435                                                 } else {
2436                                                         // Left-to-right
2437                                                         // Skip what we've matched of the overhang
2438                                                         skipLeftToRights = j;
2439                                                         if(seedOnLeft) {
2440                                                                 // ...and skip the seed if it's on the left
2441                                                                 skipLeftToRights += slen;
2442                                                         }
2443                                                 }
2444                                                 foundHit = false; // Skip this candidate
2445                                                 break;
2446                                         }
2447                                         if((int)qry[qoff] != rc) {
2448                                                 // Calculate quality of mismatched base
2449                                                 char q = phredCharToPhredQual(quals[qoff]);
2450                                                 ham += mmPenalty(this->maqPenalty_, q);
2451                                                 if(ham > this->qualMax_) {
2452                                                         // Exceeded quality limit
2453                                                         foundHit = false;
2454                                                         break;
2455                                                 }
2456                                                 // Legal mismatch outside of the anchor; record it
2457                                                 mms++;
2458                                                 nonSeedMms.push_back(qoff);
2459                                                 assert_leq(nonSeedMms.size(), (size_t)mms);
2460                                                 nonSeedRefcs.push_back("ACGTN"[rc]);
2461                                         }
2462                                 }
2463                                 if(!foundHit) continue;
2464                         }
2465                         assert(foundHit);
2466                         // Adjust ri if seed is on the right-hand side
2467                         if(!seedOnLeft) {
2468                                 ri -= readSeedOverhang;
2469                                 rir -= readSeedOverhang;
2470                         }
2471                         if(pairs != NULL) {
2472                                 TU64Pair p;
2473                                 if(ri < aoff) {
2474                                         // By convention, the upstream mate's
2475                                         // coordinates go in the 'first' field
2476                                         p.first  = ((uint64_t)tidx << 32) | (uint64_t)ri;
2477                                         p.second = ((uint64_t)tidx << 32) | (uint64_t)aoff;
2478                                 } else {
2479                                         p.second = ((uint64_t)tidx << 32) | (uint64_t)ri;
2480                                         p.first  = ((uint64_t)tidx << 32) | (uint64_t)aoff;
2481                                 }
2482                                 if(pairs->find(p) != pairs->end()) {
2483                                         // We already found this hit!  Continue.
2484                                         ASSERT_ONLY(duplicates++);
2485                                         ASSERT_ONLY(r2i++);
2486                                         continue;
2487                                 } else {
2488                                         // Record this hit
2489                                         pairs->insert(p);
2490                                 }
2491                         }
2492                         if(this->verbose_) {
2493                                 cout << "About to report seed0:" << endl;
2494                                 cout << "  ";
2495                                 for(size_t i = 0; i < qlen; i++) {
2496                                         cout << (char)qry[i];
2497                                 }
2498                                 cout << endl;
2499                                 cout << "  ";
2500                                 for(size_t i = 0; i < qlen; i++) {
2501                                         cout << "ACGT"[ref[rir+i]];
2502                                 }
2503                                 cout << endl;
2504                         }
2505                         assert_lt(r2i, r2.size());
2506                         assert_eq(re2[r2i], ri);
2507                         ranges.resize(ranges.size()+1);
2508                         Range& range = ranges.back();
2509                         assert_eq((size_t)mms, r2[r2i].numMms);
2510                         range.stratum = 0;
2511                         range.numMms = mms;
2512                         assert_eq(0, range.mms.size());
2513                         assert_eq(0, range.refcs.size());
2514                         if(mms > 0) {
2515                                 ASSERT_ONLY(size_t mmcur = 0);
2516                    &n