search_seeded_phase1.c

   1 /*
   2  * This is a fragment, included from multiple places in ebwt_search.cpp.
   3  * It implements the logic of the first 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         btf1.setReportExacts(true);
  10         bt1.setReportExacts(true);
  11
  12         if(verbose) {
  13                 cout << patFw << ":" << qual << ", " << patRc << ":" << qualRev << endl;
  14         }
  15
  16         bool done = false;
  17         if(plen < 4) {
  18                 if(!quiet) {
  19                         cerr << "Warning: Skipping read (" << name << ") because it is less than 4 characters long" << endl;
  20                 }
  21                 done = true;
  22         } else {
  23                 // Check and see if the distribution of Ns disqualifies
  24                 // this read right off the bat
  25                 size_t slen = min<size_t>(plen, seedLen);
  26                 int ns = 0;
  27                 for(size_t i = 0; i < slen; i++) {
  28                         if((int)(Dna5)patFw[i] == 4) {
  29                                 if(++ns > seedMms) {
  30                                         // Set 'done' so that
  31                                         done = true;
  32                                         break;
  33                                 }
  34                         }
  35                 }
  36         }
  37         if(done) {
  38                 ASSERT_NO_HITS_FW(true);
  39                 ASSERT_NO_HITS_RC(true);
  40                 DONEMASK_SET(patid);
  41                 skipped = true;
  42                 sink->finishRead(*patsrc, true, true);
  43                 continue;
  44         }
  45
  46         if(!nofw) {
  47                 // Do an exact-match search on the forward pattern, just in
  48                 // case we can pick it off early here
  49                 params.setFw(true);
  50                 btf1.setQuery(patsrc->bufa());
  51                 btf1.setOffs(0, plen, plen, plen, plen, plen);
  52                 if(btf1.backtrack()) {
  53                         DONEMASK_SET(patid);
  54                         continue;
  55                 }
  56         }
  57
  58         if(!norc) {
  59                 // Set up backtracker with reverse complement
  60                 params.setFw(false);
  61                 // Set up special seed bounds
  62                 if(qs < s) {
  63                         bt1.setOffs(0, 0, (seedMms > 0)? qs5 : qs,
  64                                                           (seedMms > 1)? qs5 : qs,
  65                                                           (seedMms > 2)? qs5 : qs,
  66                                                           (seedMms > 3)? qs5 : qs);
  67                 } else {
  68                         bt1.setOffs(0, 0, (seedMms > 0)? s5 : s,
  69                                                           (seedMms > 1)? s5 : s,
  70                                                           (seedMms > 2)? s5 : s,
  71                                                           (seedMms > 3)? s5 : s);
  72                 }
  73                 bt1.setQuery(patsrc->bufa());
  74                 if(bt1.backtrack()) {
  75                         // If we reach here, then we obtained a hit for case
  76                         // 1R, 2R or 3R and can stop considering this read
  77                         DONEMASK_SET(patid);
  78                         continue;
  79                 }
  80                 // If we reach here, then cases 1R, 2R, and 3R have
  81                 // been eliminated and the read needs further
  82                 // examination
  83         }
  84
  85         if(nofw && sink->finishedWithStratum(0)) { // no more exact hits are possible
  86                 DONEMASK_SET(patid);
  87                 continue;
  88         }
  89 }