2 * This is a fragment, included from multiple places in ebwt_search.cpp.
3 * It implements the logic of the second phase of the seeded, quality-
4 * aware search routine. It is implemented as a code fragment so that
5 * it can be reused in both the half-index-in-memory and full-index-in-
10 // If we reach here, then cases 1R, 2R, and 3R have been
11 // eliminated. The next most likely cases are 1F, 2F and
13 params.setFw(true); // looking at forward strand
14 btf2.setReportExacts(false);
15 btr2.setReportExacts(false);
16 btf2.setQuery(patsrc->bufa());
20 (seedMms > 0)? qs5 : qs,
21 (seedMms > 1)? qs5 : qs,
22 (seedMms > 2)? qs5 : qs,
23 (seedMms > 3)? qs5 : qs);
26 (seedMms > 0)? s5 : s,
27 (seedMms > 1)? s5 : s,
28 (seedMms > 2)? s5 : s,
29 (seedMms > 3)? s5 : s);
31 // Do a 12/24 backtrack on the forward-strand read using
32 // the mirror index. This will find all case 1F, 2F
34 if(btf2.backtrack()) {
35 // The reverse complement hit, so we're done with this
41 if(sink->finishedWithStratum(0)) { // no more exact hits are possible
47 // No need to collect partial alignments if we're not
48 // allowing mismatches in the 5' seed half
49 if(seedMms == 0) continue;
52 // If we reach here, then cases 1F, 2F, 3F, 1R, 2R, and 3R
53 // have been eliminated, leaving us with cases 4F and 4R
54 // (the cases with 1 mismatch in the 5' half of the seed)
55 params.setFw(false); // looking at reverse-comp strand
60 (seedMms > 1)? qs3 : qs,
61 (seedMms > 2)? qs3 : qs,
62 (seedMms > 3)? qs3 : qs);
66 (seedMms > 1)? s3 : s,
67 (seedMms > 2)? s3 : s,
68 (seedMms > 3)? s3 : s);
70 btr2.setQuery(patsrc->bufa());
71 btr2.setQlen(s); // just look at the seed
72 // Find partial alignments for case 4R
73 ASSERT_ONLY(bool done =) btr2.backtrack();
75 vector<PartialAlignment> partials;
76 assert(pamRc != NULL);
77 pamRc->getPartials(patid, partials);
78 if(done) assert_gt(partials.size(), 0);
79 for(size_t i = 0; i < partials.size(); i++) {
80 uint32_t pos0 = partials[i].entry.pos0;
82 uint8_t oldChar = (uint8_t)patRcRev[pos0];
83 assert_neq(oldChar, partials[i].entry.char0);
84 if(partials[i].entry.pos1 != 0xffff) {
85 uint32_t pos1 = partials[i].entry.pos1;
87 oldChar = (uint8_t)patRcRev[pos1];
88 assert_neq(oldChar, partials[i].entry.char1);
89 if(partials[i].entry.pos2 != 0xffff) {
90 uint32_t pos2 = partials[i].entry.pos2;
92 oldChar = (uint8_t)patRcRev[pos2];
93 assert_neq(oldChar, partials[i].entry.char2);