search_seeded_phase2.c

   1 /*
   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-
   6  * memory situations.
   7  */
   8 {
   9         if(!nofw) {
  10                 // If we reach here, then cases 1R, 2R, and 3R have been
  11                 // eliminated.  The next most likely cases are 1F, 2F and
  12                 // 3F...
  13                 params.setFw(true);  // looking at forward strand
  14                 btf2.setReportExacts(false);
  15                 btr2.setReportExacts(false);
  16                 btf2.setQuery(patsrc->bufa());
  17                 // Set up seed bounds
  18                 if(qs < s) {
  19                         btf2.setOffs(0, 0,
  20                                                 (seedMms > 0)? qs5 : qs,
  21                                                 (seedMms > 1)? qs5 : qs,
  22                                                 (seedMms > 2)? qs5 : qs,
  23                                                 (seedMms > 3)? qs5 : qs);
  24                 } else {
  25                         btf2.setOffs(0, 0,
  26                                                 (seedMms > 0)? s5 : s,
  27                                                 (seedMms > 1)? s5 : s,
  28                                                 (seedMms > 2)? s5 : s,
  29                                                 (seedMms > 3)? s5 : s);
  30                 }
  31                 // Do a 12/24 backtrack on the forward-strand read using
  32                 // the mirror index.  This will find all case 1F, 2F
  33                 // and 3F hits.
  34                 if(btf2.backtrack()) {
  35                         // The reverse complement hit, so we're done with this
  36                         // read
  37                         DONEMASK_SET(patid);
  38                         continue;
  39                 }
  40
  41                 if(sink->finishedWithStratum(0)) { // no more exact hits are possible
  42                         DONEMASK_SET(patid);
  43                         continue;
  44                 }
  45         }
  46
  47         // No need to collect partial alignments if we're not
  48         // allowing mismatches in the 5' seed half
  49         if(seedMms == 0) continue;
  50
  51         if(!norc) {
  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
  56                 // Set up seed bounds
  57                 if(qs < s) {
  58                         btr2.setOffs(0, 0,
  59                                                 qs3,
  60                                                 (seedMms > 1)? qs3 : qs,
  61                                                 (seedMms > 2)? qs3 : qs,
  62                                                 (seedMms > 3)? qs3 : qs);
  63                 } else {
  64                         btr2.setOffs(0, 0,
  65                                                 s3,
  66                                                 (seedMms > 1)? s3 : s,
  67                                                 (seedMms > 2)? s3 : s,
  68                                                 (seedMms > 3)? s3 : s);
  69                 }
  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();
  74 #ifndef NDEBUG
  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;
  81                         assert_lt(pos0, s5);
  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;
  86                                 assert_lt(pos1, s5);
  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;
  91                                         assert_lt(pos2, s5);
  92                                         oldChar = (uint8_t)patRcRev[pos2];
  93                                         assert_neq(oldChar, partials[i].entry.char2);
  94                                 }
  95                         }
  96                 }
  97 #endif
  98         }
  99 }