2 * This is a fragment, included from multiple places in ebwt_search.cpp.
3 * It implements the logic of the third 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-
11 btr3.setReportExacts(true);
13 btr3.setQuery(patsrc->bufa());
14 // Get all partial alignments for this read's reverse
17 if(pamRc != NULL && pamRc->size() > 0) {
18 // We can get away with an unsynchronized call because there
19 // are no writers for pamRc in this phase
20 pamRc->getPartialsUnsync(patid, pals);
22 assert_eq(0, pamRc->size());
26 // Partial alignments exist - extend them
29 btr3.setOffs(0, 0, qs, qs, qs, qs);
31 btr3.setOffs(0, 0, s, s, s, s);
33 for(size_t i = 0; i < pals.size(); i++) {
36 PartialAlignmentManager::toMutsString(
37 pals[i], patRc, qualRev, muts, !noMaqRound);
39 // Set the backtracking thresholds appropriately
40 // Now begin the backtracking, treating the first
41 // 24 bases as unrevisitable
42 ASSERT_ONLY(String<Dna5> tmp = patRc);
44 done = btr3.backtrack(oldQuals);
46 assert_eq(tmp, patRc); // assert mutations were undone
48 // The reverse complement hit, so we're done with this
51 // Got a hit; stop processing partial
55 } // Loop over partial alignments
58 // Case 4R yielded a hit continue to next pattern
60 // If we're in two-mismatch mode, then now is the time to
61 // try the final case that might apply to the reverse
62 // complement pattern: 1 mismatch in each of the 3' and 5'
63 // halves of the seed.
66 btr23.setQuery(patsrc->bufa());
67 // Set up special seed bounds
69 btr23.setOffs(qs5, qs,
71 (seedMms <= 2)? qs5 : 0, // 1revOff
72 (seedMms < 3 )? qs : qs5, // 2revOff
77 (seedMms <= 2)? s5 : 0, // 1revOff
78 (seedMms < 3 )? s : s5, // 2revOff
81 done = btr23.backtrack();
82 if(btr23.numBacktracks() == btr23.maxBacktracks()) {
87 btr23.resetNumBacktracks();
90 btr23.resetNumBacktracks();
94 if(nofw) { // no more 1-mm-in-seed hits are possible
95 //DONEMASK_SET(patid);
99 // If we reach here, then cases 1F, 2F, 3F, 1R, 2R, 3R and
100 // 4R have been eliminated leaving only 4F.
101 params.setFw(true); // looking at forward strand
102 btf3.setQuery(patsrc->bufa());
103 btf3.setQlen(seedLen); // just look at the seed
104 // Set up seed bounds
108 (seedMms > 1)? qs3 : qs,
109 (seedMms > 2)? qs3 : qs,
110 (seedMms > 3)? qs3 : qs);
114 (seedMms > 1)? s3 : s,
115 (seedMms > 2)? s3 : s,
116 (seedMms > 3)? s3 : s);
118 // Do a 12/24 seedling backtrack on the forward read
119 // using the normal index. This will find seedlings
123 vector<PartialAlignment> partials;
124 pamFw->getPartials(patid, partials);
125 for(size_t i = 0; i < partials.size(); i++) {
126 uint32_t pos0 = partials[i].entry.pos0;
128 uint8_t oldChar = (uint8_t)patFw[pos0];
129 assert_neq(oldChar, partials[i].entry.char0);
130 if(partials[i].entry.pos1 != 0xffff) {
131 uint32_t pos1 = partials[i].entry.pos1;
133 oldChar = (uint8_t)patFw[pos1];
134 assert_neq(oldChar, partials[i].entry.char1);
135 if(partials[i].entry.pos2 != 0xffff) {
136 uint32_t pos2 = partials[i].entry.pos2;
138 oldChar = (uint8_t)patFw[pos2];
139 assert_neq(oldChar, partials[i].entry.char2);