search_seeded_phase3.c

   1 /*
   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-
   6  * memory situations.
   7  */
   8 {
   9         if(!norc) {
  10                 params.setFw(false);
  11                 btr3.setReportExacts(true);
  12
  13                 btr3.setQuery(patsrc->bufa());
  14                 // Get all partial alignments for this read's reverse
  15                 // complement
  16                 pals.clear();
  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);
  21                         pamRc->clear(patid);
  22                         assert_eq(0, pamRc->size());
  23                 }
  24                 bool done = false;
  25                 if(pals.size() > 0) {
  26                         // Partial alignments exist - extend them
  27                         // Set up seed bounds
  28                         if(qs < s) {
  29                                 btr3.setOffs(0, 0, qs, qs, qs, qs);
  30                         } else {
  31                                 btr3.setOffs(0, 0, s, s, s, s);
  32                         }
  33                         for(size_t i = 0; i < pals.size(); i++) {
  34                                 seqan::clear(muts);
  35                                 uint8_t oldQuals =
  36                                         PartialAlignmentManager::toMutsString(
  37                                                 pals[i], patRc, qualRev, muts, !noMaqRound);
  38
  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);
  43                                 btr3.setMuts(&muts);
  44                                 done = btr3.backtrack(oldQuals);
  45                                 btr3.setMuts(NULL);
  46                                 assert_eq(tmp, patRc); // assert mutations were undone
  47                                 if(done) {
  48                                         // The reverse complement hit, so we're done with this
  49                                         // read
  50                                         DONEMASK_SET(patid);
  51                                         // Got a hit; stop processing partial
  52                                         // alignments
  53                                         break;
  54                                 }
  55                         } // Loop over partial alignments
  56                 }
  57                 seqan::clear(muts);
  58                 // Case 4R yielded a hit continue to next pattern
  59                 if(done) continue;
  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.
  64                 bool gaveUp = false;
  65                 if(seedMms >= 2) {
  66                         btr23.setQuery(patsrc->bufa());
  67                         // Set up special seed bounds
  68                         if(qs < s) {
  69                                 btr23.setOffs(qs5, qs,
  70                                                           0,                         // unrevOff
  71                                                           (seedMms <= 2)? qs5 : 0,   // 1revOff
  72                                                           (seedMms < 3 )? qs  : qs5, // 2revOff
  73                                                           qs);                       // 3revOff
  74                         } else {
  75                                 btr23.setOffs(s5, s,
  76                                                           0,                       // unrevOff
  77                                                           (seedMms <= 2)? s5 : 0,  // 1revOff
  78                                                           (seedMms < 3 )? s  : s5, // 2revOff
  79                                                           s);                      // 3revOff
  80                         }
  81                         done = btr23.backtrack();
  82                         if(btr23.numBacktracks() == btr23.maxBacktracks()) {
  83                                 gaveUp = true;
  84                         }
  85                         if(done) {
  86                                 DONEMASK_SET(patid);
  87                                 btr23.resetNumBacktracks();
  88                                 continue;
  89                         }
  90                         btr23.resetNumBacktracks();
  91                 }
  92         }
  93
  94         if(nofw) { // no more 1-mm-in-seed hits are possible
  95                 //DONEMASK_SET(patid);
  96                 continue;
  97         }
  98
  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
 105         if(qs < s) {
 106                 btf3.setOffs(0, 0,
 107                              qs3,
 108                              (seedMms > 1)? qs3 : qs,
 109                              (seedMms > 2)? qs3 : qs,
 110                              (seedMms > 3)? qs3 : qs);
 111         } else {
 112                 btf3.setOffs(0, 0,
 113                              s3,
 114                              (seedMms > 1)? s3 : s,
 115                              (seedMms > 2)? s3 : s,
 116                              (seedMms > 3)? s3 : s);
 117         }
 118         // Do a 12/24 seedling backtrack on the forward read
 119         // using the normal index.  This will find seedlings
 120         // for case 4F
 121         btf3.backtrack();
 122 #ifndef NDEBUG
 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;
 127                 assert_lt(pos0, s5);
 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;
 132                         assert_lt(pos1, s5);
 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;
 137                                 assert_lt(pos2, s5);
 138                                 oldChar = (uint8_t)patFw[pos2];
 139                                 assert_neq(oldChar, partials[i].entry.char2);
 140                         }
 141                 }
 142         }
 143 #endif
 144 }