aligner.h

   1 /**
   2  * aligner.h
   3  *
   4  * A generic class providing a stateful way to find alignments.
   5  */
   6
   7 #ifndef ALIGNER_H_
   8 #define ALIGNER_H_
   9
  10 #include <iostream>
  11 #include <set>
  12 #include <stdint.h>
  13 #include "seqan/sequence.h"
  14 #include "assert_helpers.h"
  15 #include "ebwt.h"
  16 #include "pat.h"
  17 #include "range.h"
  18 #include "range_source.h"
  19 #include "range_chaser.h"
  20 #include "ref_aligner.h"
  21 #include "reference.h"
  22 #include "aligner_metrics.h"
  23 #include "search_globals.h"
  24
  25 /**
  26  * State machine for carrying out an alignment, which usually consists
  27  * of a series of phases that conduct different alignments using
  28  * different backtracking constraints.
  29  *
  30  * Each Aligner should have a dedicated PatternSourcePerThread.
  31  */
  32 class Aligner {
  33 public:
  34         Aligner(bool _done, bool rangeMode) :
  35                 done(_done), patsrc_(NULL), bufa_(NULL), bufb_(NULL),
  36                 rangeMode_(rangeMode)
  37         { }
  38
  39         virtual ~Aligner() { }
  40         /// Advance the range search by one memory op
  41         virtual bool advance() = 0;
  42
  43         /// Prepare Aligner for the next read
  44         virtual void setQuery(PatternSourcePerThread *patsrc) {
  45                 assert(patsrc != NULL);
  46                 patsrc_ = patsrc;
  47                 bufa_ = &patsrc->bufa();
  48                 assert(bufa_ != NULL);
  49                 bufb_ = &patsrc->bufb();
  50                 alen_ = bufa_->length();
  51                 blen_ = (bufb_ != NULL) ? bufb_->length() : 0;
  52                 rand_.init(bufa_->seed);
  53         }
  54
  55         /**
  56          * Set to true if all searching w/r/t the current query is
  57          * finished or if there is no current query.
  58          */
  59         bool done;
  60
  61 protected:
  62
  63         // Current read pair
  64         PatternSourcePerThread* patsrc_;
  65         ReadBuf* bufa_;
  66         uint32_t alen_;
  67         ReadBuf* bufb_;
  68         uint32_t blen_;
  69         bool rangeMode_;
  70         RandomSource rand_;
  71 };
  72
  73 /**
  74  * Abstract parent factory class for constructing aligners of all kinds.
  75  */
  76 class AlignerFactory {
  77 public:
  78         virtual ~AlignerFactory() { }
  79         virtual Aligner* create() const = 0;
  80
  81         /**
  82          * Allocate a vector of n Aligners; use destroy(std::vector...) to
  83          * free the memory.
  84          */
  85         virtual std::vector<Aligner*>* create(uint32_t n) const {
  86                 std::vector<Aligner*>* v = new std::vector<Aligner*>;
  87                 for(uint32_t i = 0; i < n; i++) {
  88                         v->push_back(create());
  89                         assert(v->back() != NULL);
  90                 }
  91                 return v;
  92         }
  93
  94         /// Free memory associated with the aligner
  95         virtual void destroy(Aligner* al) const {
  96                 assert(al != NULL);
  97                 // Free the Aligner
  98                 delete al;
  99         }
 100
 101         /// Free memory associated with an aligner list
 102         virtual void destroy(std::vector<Aligner*>* als) const {
 103                 assert(als != NULL);
 104                 // Free all of the Aligners
 105                 for(size_t i = 0; i < als->size(); i++) {
 106                         if((*als)[i] != NULL) {
 107                                 delete (*als)[i];
 108                                 (*als)[i] = NULL;
 109                         }
 110                 }
 111                 // Free the vector
 112                 delete als;
 113         }
 114 };
 115
 116 /**
 117  * Coordinates multiple aligners of the same type (i.e. either all
 118  * single-end or all paired-end).
 119  */
 120 class MultiAligner {
 121 public:
 122         MultiAligner(
 123                         uint32_t n,
 124                         uint32_t qUpto,
 125                         const AlignerFactory& alignFact,
 126                         const PatternSourcePerThreadFactory& patsrcFact) :
 127                         n_(n), qUpto_(qUpto),
 128                         alignFact_(alignFact), patsrcFact_(patsrcFact),
 129                         aligners_(NULL), patsrcs_(NULL)
 130         {
 131                 aligners_ = alignFact_.create(n_);
 132                 assert(aligners_ != NULL);
 133                 patsrcs_ = patsrcFact_.create(n_);
 134                 assert(patsrcs_ != NULL);
 135         }
 136
 137         /// Free memory associated with the aligners and their pattern sources.
 138         virtual ~MultiAligner() {
 139                 alignFact_.destroy(aligners_);
 140                 patsrcFact_.destroy(patsrcs_);
 141         }
 142
 143         /**
 144          * Advance an array of aligners in parallel, using prefetches to
 145          * try to hide all the latency.
 146          */
 147         void run() {
 148                 bool done = false;
 149                 while(!done) {
 150                         done = true;
 151                         for(uint32_t i = 0; i < n_; i++) {
 152                                 if(!(*aligners_)[i]->done) {
 153                                         // Advance an aligner already in progress
 154                                         done = false;
 155                                         (*aligners_)[i]->advance();
 156                                 } else {
 157                                         // Get a new read and initialize an aligner with it
 158                                         (*patsrcs_)[i]->nextReadPair();
 159                                         if(!(*patsrcs_)[i]->empty() && (*patsrcs_)[i]->patid() < qUpto_) {
 160                                                 (*aligners_)[i]->setQuery((*patsrcs_)[i]);
 161                                                 assert(!(*aligners_)[i]->done);
 162                                                 done = false;
 163                                         } else {
 164                                                 // No more reads; if done == true, it remains
 165                                                 // true
 166                                         }
 167                                 }
 168                         }
 169                 }
 170         }
 171
 172 protected:
 173         uint32_t n_;     /// Number of aligners
 174         uint32_t qUpto_; /// Number of reads to align before stopping
 175         const AlignerFactory&                  alignFact_;
 176         const PatternSourcePerThreadFactory&   patsrcFact_;
 177         std::vector<Aligner *>*                aligners_;
 178         std::vector<PatternSourcePerThread *>* patsrcs_;
 179 };
 180
 181 /**
 182  * Coordinates multiple single-end and paired-end aligners, routing
 183  * reads to one or the other type as appropriate.
 184  */
 185 class MixedMultiAligner {
 186 public:
 187         MixedMultiAligner(
 188                         uint32_t n,
 189                         uint32_t qUpto,
 190                         const AlignerFactory& alignSEFact,
 191                         const AlignerFactory& alignPEFact,
 192                         const PatternSourcePerThreadFactory& patsrcFact) :
 193                         n_(n), qUpto_(qUpto),
 194                         alignSEFact_(alignSEFact),
 195                         alignPEFact_(alignPEFact),
 196                         patsrcFact_(patsrcFact),
 197                         alignersSE_(NULL),
 198                         alignersPE_(NULL),
 199                         seOrPe_(NULL),
 200                         patsrcs_(NULL)
 201         {
 202                 // Instantiate all single-end aligners
 203                 alignersSE_ = alignSEFact_.create(n_);
 204                 assert(alignersSE_ != NULL);
 205                 // Instantiate all paired-end aligners
 206                 alignersPE_ = alignPEFact_.create(n_);
 207                 assert(alignersPE_ != NULL);
 208                 // Allocate array of boolean flags indicating whether each of
 209                 // the slots is currently using the single-end or paired-end
 210                 // aligner
 211                 seOrPe_ = new bool[n_];
 212                 for(uint32_t i = 0; i < n_; i++) {
 213                         seOrPe_[i] = true;
 214                 }
 215                 // Instantiate all read sources
 216                 patsrcs_ = patsrcFact_.create(n_);
 217                 assert(patsrcs_ != NULL);
 218         }
 219
 220         /// Free memory associated with the aligners and their pattern sources.
 221         virtual ~MixedMultiAligner() {
 222                 alignSEFact_.destroy(alignersSE_);
 223                 alignPEFact_.destroy(alignersPE_);
 224                 patsrcFact_.destroy(patsrcs_);
 225                 delete[] seOrPe_;
 226         }
 227
 228         /**
 229          * Advance an array of aligners in parallel, using prefetches to
 230          * try to hide all the latency.
 231          */
 232         void run(bool verbose = false) {
 233                 bool done = false;
 234                 bool first = true;
 235                 if(n_ == 1) {
 236                         Aligner *al = seOrPe_[0] ? (*alignersSE_)[0] : (*alignersPE_)[0];
 237                         PatternSourcePerThread *ps = (*patsrcs_)[0];
 238                         while(!done) {
 239                                 done = true;
 240                                 if(!first && !al->done) {
 241                                         // Advance an aligner already in progress; this is
 242                                         // the common case
 243                                         done = false;
 244                                         al->advance();
 245                                 } else {
 246                                         // Get a new read
 247                                         ps->nextReadPair();
 248                                         if(ps->patid() < qUpto_ && !ps->empty()) {
 249                                                 if(ps->paired()) {
 250                                                         // Read currently in buffer is paired-end
 251                                                         (*alignersPE_)[0]->setQuery(ps);
 252                                                         al = (*alignersPE_)[0];
 253                                                         seOrPe_[0] = false; // false -> paired
 254                                                 } else {
 255                                                         // Read currently in buffer is single-end
 256                                                         (*alignersSE_)[0]->setQuery(ps);
 257                                                         al = (*alignersSE_)[0];
 258                                                         seOrPe_[0] = true; // true = unpaired
 259                                                 }
 260                                                 done = false;
 261                                         } else {
 262                                                 // No more reads; if done == true, it remains
 263                                                 // true
 264                                         }
 265                                 }
 266                                 first = false;
 267                         }
 268                 } else {
 269                         while(!done) {
 270                                 done = true;
 271                                 for(uint32_t i = 0; i < n_; i++) {
 272                                         Aligner *al = seOrPe_[i] ? (*alignersSE_)[i] :
 273                                                                                            (*alignersPE_)[i];
 274                                         if(!first && !al->done) {
 275                                                 // Advance an aligner already in progress; this is
 276                                                 // the common case
 277                                                 done = false;
 278                                                 al->advance();
 279                                         } else {
 280                                                 // Feed a new read to a vacant aligner
 281                                                 PatternSourcePerThread *ps = (*patsrcs_)[i];
 282                                                 // Get a new read
 283                                                 ps->nextReadPair();
 284                                                 if(ps->patid() < qUpto_ && !ps->empty()) {
 285                                                         if(ps->paired()) {
 286                                                                 // Read currently in buffer is paired-end
 287                                                                 (*alignersPE_)[i]->setQuery(ps);
 288                                                                 seOrPe_[i] = false; // false -> paired
 289                                                         } else {
 290                                                                 // Read currently in buffer is single-end
 291                                                                 (*alignersSE_)[i]->setQuery(ps);
 292                                                                 seOrPe_[i] = true; // true = unpaired
 293                                                         }
 294                                                         done = false;
 295                                                 } else {
 296                                                         // No more reads; if done == true, it remains
 297                                                         // true
 298                                                 }
 299                                         }
 300                                 }
 301                                 first = false;
 302                         }
 303                 }
 304         }
 305
 306 protected:
 307         uint32_t n_;     /// Number of aligners
 308         uint32_t qUpto_; /// Number of reads to align before stopping
 309         const AlignerFactory&                  alignSEFact_;
 310         const AlignerFactory&                  alignPEFact_;
 311         const PatternSourcePerThreadFactory&   patsrcFact_;
 312         std::vector<Aligner *>*                alignersSE_;
 313         std::vector<Aligner *>*                alignersPE_;
 314         bool *                                 seOrPe_;
 315         std::vector<PatternSourcePerThread *>* patsrcs_;
 316 };
 317
 318 /**
 319  * An aligner for finding exact matches of unpaired reads.  Always
 320  * tries the forward-oriented version of the read before the reverse-
 321  * oriented read.
 322  */
 323 template<typename TRangeSource>
 324 class UnpairedAlignerV2 : public Aligner {
 325         typedef RangeSourceDriver<TRangeSource> TDriver;
 326 public:
 327         UnpairedAlignerV2(
 328                 EbwtSearchParams<String<Dna> >* params,
 329                 TDriver* driver,
 330                 RangeChaser<String<Dna> >* rchase,
 331                 HitSink& sink,
 332                 const HitSinkPerThreadFactory& sinkPtFactory,
 333                 HitSinkPerThread* sinkPt,
 334                 vector<String<Dna5> >& os,
 335                 const BitPairReference* refs,
 336                 bool rangeMode,
 337                 bool verbose,
 338                 bool quiet,
 339                 int maxBts,
 340                 ChunkPool *pool,
 341                 int *btCnt = NULL,
 342                 AlignerMetrics *metrics = NULL) :
 343                 Aligner(true, rangeMode),
 344                 refs_(refs),
 345                 doneFirst_(true),
 346                 firstIsFw_(true),
 347                 chase_(false),
 348                 sinkPtFactory_(sinkPtFactory),
 349                 sinkPt_(sinkPt),
 350                 params_(params),
 351                 rchase_(rchase),
 352                 driver_(driver),
 353                 verbose_(verbose),
 354                 quiet_(quiet),
 355                 maxBts_(maxBts),
 356                 pool_(pool),
 357                 btCnt_(btCnt),
 358                 metrics_(metrics)
 359         {
 360                 assert(pool_   != NULL);
 361                 assert(sinkPt_ != NULL);
 362                 assert(params_ != NULL);
 363                 assert(driver_ != NULL);
 364         }
 365
 366         virtual ~UnpairedAlignerV2() {
 367                 delete driver_;  driver_  = NULL;
 368                 delete params_;  params_  = NULL;
 369                 delete rchase_;  rchase_  = NULL;
 370                 delete[] btCnt_; btCnt_   = NULL;
 371                 sinkPtFactory_.destroy(sinkPt_); sinkPt_ = NULL;
 372         }
 373
 374         /**
 375          * Prepare this aligner for the next read.
 376          */
 377         virtual void setQuery(PatternSourcePerThread* patsrc) {
 378                 Aligner::setQuery(patsrc); // set fields & random seed
 379                 if(metrics_ != NULL) {
 380                         metrics_->nextRead(patsrc->bufa().patFw);
 381                 }
 382                 pool_->reset(&patsrc->bufa().name, patsrc->patid());
 383                 if(patsrc->bufa().length() < 4) {
 384                         if(!quiet_) {
 385                                 cerr << "Warning: Skipping read " << patsrc->bufa().name
 386                                      << " because it is less than 4 characters long" << endl;
 387                         }
 388                         this->done = true;
 389                         sinkPt_->finishRead(*patsrc_, true, true);
 390                         return;
 391                 }
 392                 driver_->setQuery(patsrc, NULL);
 393                 this->done = driver_->done;
 394                 doneFirst_ = false;
 395                 // Reset #-backtrack countdown
 396                 if(btCnt_ != NULL) *btCnt_ = maxBts_;
 397                 if(sinkPt_->setHits(patsrc->bufa().hitset)) {
 398                         this->done = true;
 399                         sinkPt_->finishRead(*patsrc_, true, true);
 400                 }
 401                 // Grab a bit from the pseudo-random seed to determine whether
 402                 // to start with forward or reverse complement
 403                 firstIsFw_ = ((patsrc->bufa().seed & 0x10) == 0);
 404                 chase_ = false;
 405         }
 406
 407         /**
 408          * Helper for reporting an alignment.
 409          */
 410         inline bool report(const Range& ra,
 411                            uint32_t first,
 412                            uint32_t second,
 413                            uint32_t tlen)
 414         {
 415                 bool ebwtFw = ra.ebwt->fw();
 416                 params_->setFw(ra.fw);
 417                 assert_eq(bufa_->color, color);
 418                 return params_->reportHit(
 419                                 ra.fw ? (ebwtFw? bufa_->patFw    : bufa_->patFwRev) :
 420                                         (ebwtFw? bufa_->patRc    : bufa_->patRcRev),
 421                                 ra.fw ? (ebwtFw? &bufa_->qual    : &bufa_->qualRev) :
 422                                         (ebwtFw? &bufa_->qualRev : &bufa_->qual),
 423                                 &bufa_->name,
 424                                 bufa_->color,
 425                                 colorExEnds,
 426                                 snpPhred,
 427                                 refs_,
 428                                 ra.ebwt->rmap(),
 429                                 ebwtFw,
 430                                 ra.mms,                   // mismatch positions
 431                                 ra.refcs,                 // reference characters for mms
 432                                 ra.numMms,                // # mismatches
 433                                 make_pair(first, second), // position
 434                                 make_pair(0, 0),          // (bogus) mate position
 435                                 true,                     // (bogus) mate orientation
 436                                 0,                        // (bogus) mate length
 437                                 make_pair(ra.top, ra.bot),// arrows
 438                                 tlen,                     // textlen
 439                                 alen_,                    // qlen
 440                                 ra.stratum,               // alignment stratum
 441                                 ra.cost,                  // cost, including qual penalty
 442                                 ra.bot - ra.top - 1,      // # other hits
 443                                 patsrc_->patid(),         // pattern id
 444                                 bufa_->seed,              // pseudo-random seed
 445                                 0);                       // mate (0 = unpaired)
 446         }
 447
 448         /**
 449          * Advance the aligner.  Return true iff we're
 450          * done with this read.
 451          */
 452         virtual bool advance() {
 453                 assert(!this->done);
 454                 if(chase_) {
 455                         assert(!rangeMode_);
 456                         assert(driver_->foundRange);
 457                         assert(!sinkPt_->irrelevantCost(driver_->range().cost));
 458                         if(!rchase_->foundOff() && !rchase_->done) {
 459                                 rchase_->advance();
 460                                 return false;
 461                         }
 462                         if(rchase_->foundOff()) {
 463                                 this->done = report(driver_->range(), rchase_->off().first,
 464                                                     rchase_->off().second, rchase_->tlen());
 465                                 rchase_->reset();
 466                         } else {
 467                                 assert(rchase_->done);
 468                                 // Forget this range; keep looking for ranges
 469                                 chase_ = false;
 470                                 driver_->foundRange = false;
 471                                 this->done = driver_->done;
 472                         }
 473                 }
 474                 // Still advancing a
 475                 if(!this->done && !chase_) {
 476                         assert(!driver_->done || driver_->foundRange);
 477                         if(driver_->foundRange) {
 478                                 const Range& ra = driver_->range();
 479                                 assert(!sinkPt_->irrelevantCost(ra.cost));
 480                                 assert(ra.repOk());
 481                                 if(rangeMode_) {
 482                                         this->done = report(ra, ra.top, ra.bot, 0);
 483                                         driver_->foundRange = false;
 484                                 } else {
 485                                         rchase_->setTopBot(ra.top, ra.bot, alen_, rand_, ra.ebwt);
 486                                         if(rchase_->foundOff()) {
 487                                                 this->done = report(
 488                                                                 ra, rchase_->off().first,
 489                                                                 rchase_->off().second, rchase_->tlen());
 490                                                 rchase_->reset();
 491                                         }
 492                                         if(!rchase_->done && !sinkPt_->irrelevantCost(ra.cost)) {
 493                                                 // Keep chasing this range
 494                                                 chase_ = true;
 495                                         } else {
 496                                                 driver_->foundRange = false;
 497                                         }
 498                                 }
 499                         } else {
 500                                 this->done = sinkPt_->irrelevantCost(driver_->minCost);
 501                                 if(!this->done) {
 502                                         driver_->advance(ADV_COST_CHANGES);
 503                                 } else {
 504                                         // No longer necessarily true with chain input
 505                                         //assert(!sinkPt_->spanStrata());
 506                                 }
 507                         }
 508                         if(driver_->done && !driver_->foundRange && !chase_) {
 509                                 this->done = true;
 510                         }
 511                 }
 512                 if(this->done) {
 513                         sinkPt_->finishRead(*patsrc_, true, true);
 514                 }
 515                 return this->done;
 516         }
 517
 518 protected:
 519
 520         // Reference sequences (needed for colorspace decoding)
 521         const BitPairReference* refs_;
 522
 523         // Progress state
 524         bool doneFirst_;
 525         bool firstIsFw_;
 526         bool chase_;
 527
 528         // Temporary HitSink; to be deleted
 529         const HitSinkPerThreadFactory& sinkPtFactory_;
 530         HitSinkPerThread* sinkPt_;
 531
 532         // State for alignment
 533         EbwtSearchParams<String<Dna> >* params_;
 534
 535         // State for getting alignments from ranges statefully
 536         RangeChaser<String<Dna> >* rchase_;
 537
 538         // Range-finding state
 539         TDriver* driver_;
 540
 541         bool verbose_; // be talkative
 542         bool quiet_; // don't print informational/warning info
 543
 544         const int maxBts_;
 545         ChunkPool *pool_;
 546         int *btCnt_;
 547         AlignerMetrics *metrics_;
 548 };
 549
 550 /**
 551  * An aligner for finding paired alignments while operating entirely
 552  * within the Burrows-Wheeler domain.
 553  */
 554 template<typename TRangeSource>
 555 class PairedBWAlignerV1 : public Aligner {
 556
 557         typedef std::pair<uint32_t,uint32_t> U32Pair;
 558         typedef std::vector<U32Pair> U32PairVec;
 559         typedef std::vector<Range> TRangeVec;
 560         typedef RangeSourceDriver<TRangeSource> TDriver;
 561         typedef std::pair<uint64_t, uint64_t> TU64Pair;
 562         typedef std::set<TU64Pair> TSetPairs;
 563
 564 public:
 565         PairedBWAlignerV1(
 566                 EbwtSearchParams<String<Dna> >* params,
 567                 TDriver* driver1Fw, TDriver* driver1Rc,
 568                 TDriver* driver2Fw, TDriver* driver2Rc,
 569                 RefAligner<String<Dna5> >* refAligner,
 570                 RangeChaser<String<Dna> >* rchase,
 571                 HitSink& sink,
 572                 const HitSinkPerThreadFactory& sinkPtFactory,
 573                 HitSinkPerThread* sinkPt,
 574                 bool fw1, bool fw2,
 575                 uint32_t minInsert,
 576                 uint32_t maxInsert,
 577                 bool dontReconcile,
 578                 uint32_t symCeiling,
 579                 uint32_t mixedThresh,
 580                 uint32_t mixedAttemptLim,
 581                 const BitPairReference* refs,
 582                 bool rangeMode,
 583                 bool verbose,
 584                 bool quiet,
 585                 int maxBts,
 586                 ChunkPool *pool,
 587                 int *btCnt) :
 588                 Aligner(true, rangeMode),
 589                 refs_(refs),
 590                 patsrc_(NULL), qlen1_(0), qlen2_(0), doneFw_(true),
 591                 doneFwFirst_(true),
 592                 chase1Fw_(false), chase1Rc_(false),
 593                 chase2Fw_(false), chase2Rc_(false),
 594                 delayedChase1Fw_(false), delayedChase1Rc_(false),
 595                 delayedChase2Fw_(false), delayedChase2Rc_(false),
 596                 refAligner_(refAligner),
 597                 sinkPtFactory_(sinkPtFactory),
 598                 sinkPt_(sinkPt),
 599                 params_(params),
 600                 minInsert_(minInsert),
 601                 maxInsert_(maxInsert),
 602                 dontReconcile_(dontReconcile),
 603                 symCeiling_(symCeiling),
 604                 mixedThresh_(mixedThresh),
 605                 mixedAttemptLim_(mixedAttemptLim),
 606                 mixedAttempts_(0),
 607                 fw1_(fw1), fw2_(fw2),
 608                 rchase_(rchase),
 609                 verbose_(verbose),
 610                 quiet_(quiet),
 611                 maxBts_(maxBts),
 612                 pool_(pool),
 613                 btCnt_(btCnt),
 614                 driver1Fw_(driver1Fw), driver1Rc_(driver1Rc),
 615                 offs1FwSz_(0), offs1RcSz_(0),
 616                 driver2Fw_(driver2Fw), driver2Rc_(driver2Rc),
 617                 offs2FwSz_(0), offs2RcSz_(0),
 618
 619                 chaseL_fw_       (fw1_ ? chase1Fw_        : chase1Rc_),
 620                 chaseR_fw_       (fw2_ ? chase2Fw_        : chase2Rc_),
 621                 delayedchaseL_fw_(fw1_ ? delayedChase1Fw_ : delayedChase1Rc_),
 622                 delayedchaseR_fw_(fw2_ ? delayedChase2Fw_ : delayedChase2Rc_),
 623                 drL_fw_          (fw1_ ? *driver1Fw_      : *driver1Rc_),
 624                 drR_fw_          (fw2_ ? *driver2Fw_      : *driver2Rc_),
 625                 offsLarr_fw_     (fw1_ ? offs1FwArr_      : offs1RcArr_),
 626                 offsRarr_fw_     (fw2_ ? offs2FwArr_      : offs2RcArr_),
 627                 rangesLarr_fw_   (fw1_ ? ranges1FwArr_    : ranges1RcArr_),
 628                 rangesRarr_fw_   (fw2_ ? ranges2FwArr_    : ranges2RcArr_),
 629                 offsLsz_fw_      (fw1_ ? offs1FwSz_       : offs1RcSz_),
 630                 offsRsz_fw_      (fw2_ ? offs2FwSz_       : offs2RcSz_),
 631
 632                 chaseL_rc_       (fw2_ ? chase2Rc_        : chase2Fw_),
 633                 chaseR_rc_       (fw1_ ? chase1Rc_        : chase1Fw_),
 634                 delayedchaseL_rc_(fw2_ ? delayedChase2Rc_ : delayedChase2Fw_),
 635                 delayedchaseR_rc_(fw1_ ? delayedChase1Rc_ : delayedChase1Fw_),
 636                 drL_rc_          (fw2_ ? *driver2Rc_      : *driver2Fw_),
 637                 drR_rc_          (fw1_ ? *driver1Rc_      : *driver1Fw_),
 638                 offsLarr_rc_     (fw2_ ? offs2RcArr_      : offs2FwArr_),
 639                 offsRarr_rc_     (fw1_ ? offs1RcArr_      : offs1FwArr_),
 640                 rangesLarr_rc_   (fw2_ ? ranges2RcArr_    : ranges2FwArr_),
 641                 rangesRarr_rc_   (fw1_ ? ranges1RcArr_    : ranges1FwArr_),
 642                 offsLsz_rc_      (fw2_ ? offs2RcSz_       : offs2FwSz_),
 643                 offsRsz_rc_      (fw1_ ? offs1RcSz_       : offs1FwSz_),
 644
 645                 chaseL_       (&chaseL_fw_),
 646                 chaseR_       (&chaseR_fw_),
 647                 delayedchaseL_(&delayedchaseL_fw_),
 648                 delayedchaseR_(&delayedchaseR_fw_),
 649                 drL_          (&drL_fw_),
 650                 drR_          (&drR_fw_),
 651                 offsLarr_     (offsLarr_fw_),
 652                 offsRarr_     (offsRarr_fw_),
 653                 rangesLarr_   (rangesLarr_fw_),
 654                 rangesRarr_   (rangesRarr_fw_),
 655                 offsLsz_      (&offsLsz_fw_),
 656                 offsRsz_      (&offsRsz_fw_),
 657                 donePair_     (&doneFw_),
 658                 fwL_(fw1),
 659                 fwR_(fw2),
 660                 verbose2_(false)
 661         {
 662                 assert(pool_      != NULL);
 663                 assert(sinkPt_    != NULL);
 664                 assert(params_    != NULL);
 665                 assert(driver1Fw_ != NULL);
 666                 assert(driver1Rc_ != NULL);
 667                 assert(driver2Fw_ != NULL);
 668                 assert(driver2Rc_ != NULL);
 669         }
 670
 671         virtual ~PairedBWAlignerV1() {
 672                 delete driver1Fw_; driver1Fw_ = NULL;
 673                 delete driver1Rc_; driver1Rc_ = NULL;
 674                 delete driver2Fw_; driver2Fw_ = NULL;
 675                 delete driver2Rc_; driver2Rc_ = NULL;
 676                 delete params_;    params_    = NULL;
 677                 delete rchase_;    rchase_    = NULL;
 678                 delete[] btCnt_;   btCnt_     = NULL;
 679                 delete refAligner_; refAligner_ = NULL;
 680                 sinkPtFactory_.destroy(sinkPt_); sinkPt_ = NULL;
 681         }
 682
 683         /**
 684          * Prepare this aligner for the next read.
 685          */
 686         virtual void setQuery(PatternSourcePerThread* patsrc) {
 687                 assert(!patsrc->bufa().empty());
 688                 Aligner::setQuery(patsrc); // set fields & random seed
 689                 assert(!patsrc->bufb().empty());
 690                 // Give all of the drivers pointers to the relevant read info
 691                 patsrc_ = patsrc;
 692                 pool_->reset(&patsrc->bufa().name, patsrc->patid());
 693                 if(patsrc->bufa().length() < 4 || patsrc->bufb().length() < 4) {
 694                         if(!quiet_) {
 695                                 cerr << "Warning: Skipping pair " << patsrc->bufa().name
 696                                          << " because a mate is less than 4 characters long" << endl;
 697                         }
 698                         this->done = true;
 699                         sinkPt_->finishRead(*patsrc_, true, true);
 700                         return;
 701                 }
 702                 driver1Fw_->setQuery(patsrc, NULL);
 703                 driver1Rc_->setQuery(patsrc, NULL);
 704                 driver2Fw_->setQuery(patsrc, NULL);
 705                 driver2Rc_->setQuery(patsrc, NULL);
 706                 qlen1_ = patsrc_->bufa().length();
 707                 qlen2_ = patsrc_->bufb().length();
 708                 if(btCnt_ != NULL) (*btCnt_) = maxBts_;
 709                 // Neither orientation is done
 710                 doneFw_   = false;
 711                 doneFwFirst_ = true;
 712                 this->done   = false;
 713                 // No ranges are being chased yet
 714                 chase1Fw_ = false;
 715                 chase1Rc_ = false;
 716                 chase2Fw_ = false;
 717                 chase2Rc_ = false;
 718                 delayedChase1Fw_ = false;
 719                 delayedChase1Rc_ = false;
 720                 delayedChase2Fw_ = false;
 721                 delayedChase2Rc_ = false;
 722                 // Clear all intermediate ranges
 723                 for(size_t i = 0; i < 32; i++) {
 724                         offs1FwArr_[i].clear();   offs1RcArr_[i].clear();
 725                         offs2FwArr_[i].clear();   offs2RcArr_[i].clear();
 726                         ranges1FwArr_[i].clear(); ranges1RcArr_[i].clear();
 727                         ranges2FwArr_[i].clear(); ranges2RcArr_[i].clear();
 728                 }
 729                 offs1FwSz_ = offs1RcSz_ = offs2FwSz_ = offs2RcSz_ = 0;
 730                 chaseL_        = &chaseL_fw_;
 731                 chaseR_        = &chaseR_fw_;
 732                 delayedchaseL_ = &delayedchaseL_fw_;
 733                 delayedchaseR_ = &delayedchaseR_fw_;
 734                 drL_           = &drL_fw_;
 735                 drR_           = &drR_fw_;
 736                 offsLarr_      = offsLarr_fw_;
 737                 offsRarr_      = offsRarr_fw_;
 738                 rangesLarr_    = rangesLarr_fw_;
 739                 rangesRarr_    = rangesRarr_fw_;
 740                 offsLsz_       = &offsLsz_fw_;
 741                 offsRsz_       = &offsRsz_fw_;
 742                 donePair_      = &doneFw_;
 743                 fwL_           = fw1_;
 744                 fwR_           = fw2_;
 745                 mixedAttempts_ = 0;
 746                 pairs_fw_.clear();
 747                 pairs_rc_.clear();
 748 #ifndef NDEBUG
 749                 allTopsL_fw_.clear();
 750                 allTopsR_fw_.clear();
 751                 allTopsL_rc_.clear();
 752                 allTopsR_rc_.clear();
 753 #endif
 754         }
 755
 756         /**
 757          * Advance the aligner by one memory op.  Return true iff we're
 758          * done with this read.
 759          *
 760          * A call to this function does one of many things:
 761          * 1. Advance a RangeSourceDriver and check if it found a new range
 762          * 2. Advance a RowChaseDriver and check if it found a reference
 763          *    offset for a an alignment in a range
 764          */
 765         virtual bool advance() {
 766                 assert(!this->done);
 767                 if(doneFw_ && doneFwFirst_) {
 768                         if(verbose2_) cout << "--" << endl;
 769                         chaseL_        = &chaseL_rc_;
 770                         chaseR_        = &chaseR_rc_;
 771                         delayedchaseL_ = &delayedchaseL_rc_;
 772                         delayedchaseR_ = &delayedchaseR_rc_;
 773                         drL_           = &drL_rc_;
 774                         drR_           = &drR_rc_;
 775                         offsLarr_      = offsLarr_rc_;
 776                         offsRarr_      = offsRarr_rc_;
 777                         rangesLarr_    = rangesLarr_rc_;
 778                         rangesRarr_    = rangesRarr_rc_;
 779                         offsLsz_       = &offsLsz_rc_;
 780                         offsRsz_       = &offsRsz_rc_;
 781                         donePair_      = &this->done;
 782                         fwL_           = !fw2_;
 783                         fwR_           = !fw1_;
 784                         doneFwFirst_   = false;
 785                         mixedAttempts_ = 0;
 786                 }
 787                 bool chasing = *chaseL_ || *chaseR_;
 788                 if(chasing && !rchase_->foundOff() && !rchase_->done) {
 789                         rchase_->advance();
 790                         return false;
 791                 }
 792                 advanceOrientation(!doneFw_);
 793                 if(this->done) {
 794                         if(verbose2_) cout << "----" << endl;
 795                         sinkPt_->finishRead(*patsrc_, true, true);
 796                 }
 797                 return this->done;
 798         }
 799
 800 protected:
 801
 802         /**
 803          * Helper for reporting a pair of alignments.  As of now, we report
 804          * a paired alignment by reporting two consecutive alignments, one
 805          * for each mate.
 806          */
 807         bool report(const Range& rL, // range for upstream mate
 808                     const Range& rR, // range for downstream mate
 809                     uint32_t first,  // ref idx
 810                     uint32_t upstreamOff, // offset for upstream mate
 811                     uint32_t dnstreamOff, // offset for downstream mate
 812                     uint32_t tlen, // length of ref
 813                     bool pairFw,   // whether the pair is being mapped to fw strand
 814                     bool ebwtFwL,
 815                     bool ebwtFwR,
 816                     const ReferenceMap* rmap)
 817         {
 818                 assert_lt(upstreamOff, dnstreamOff);
 819                 uint32_t spreadL = rL.bot - rL.top;
 820                 uint32_t spreadR = rR.bot - rR.top;
 821                 uint32_t oms = min(spreadL, spreadR) - 1;
 822                 ReadBuf* bufL = pairFw ? bufa_ : bufb_;
 823                 ReadBuf* bufR = pairFw ? bufb_ : bufa_;
 824                 uint32_t lenL = pairFw ? alen_ : blen_;
 825                 uint32_t lenR = pairFw ? blen_ : alen_;
 826                 bool ret;
 827                 assert(!params_->sink().exceededOverThresh());
 828                 params_->setFw(rL.fw);
 829                 assert_eq(bufL->color, color);
 830                 // Print upstream mate first
 831                 ret = params_->reportHit(
 832                                 rL.fw ? (ebwtFwL?  bufL->patFw  :  bufL->patFwRev) :
 833                                             (ebwtFwL?  bufL->patRc  :  bufL->patRcRev),
 834                                 rL.fw ? (ebwtFwL? &bufL->qual    : &bufL->qualRev) :
 835                                         (ebwtFwL? &bufL->qualRev : &bufL->qual),
 836                                 &bufL->name,
 837                                 bufL->color,
 838                                 colorExEnds,
 839                                 snpPhred,
 840                                 refs_,
 841                                 rmap,
 842                                 ebwtFwL,
 843                                 rL.mms,                       // mismatch positions
 844                                 rL.refcs,                     // reference characters for mms
 845                                 rL.numMms,                    // # mismatches
 846                                 make_pair(first, upstreamOff),// position
 847                                 make_pair(first, dnstreamOff),// mate position
 848                                 rR.fw,                        // mate orientation
 849                                 lenR,                         // mate length
 850                                 make_pair(rL.top, rL.bot),    // arrows
 851                                 tlen,                         // textlen
 852                                 lenL,                         // qlen
 853                                 rL.stratum,                   // alignment stratum
 854                                 rL.cost,                      // cost, including quality penalty
 855                                 oms,                          // # other hits
 856                                 bufL->patid,
 857                                 bufL->seed,
 858                                 pairFw ? 1 : 2);
 859                 if(ret) {
 860                         return true; // can happen when -m is set
 861                 }
 862                 params_->setFw(rR.fw);
 863                 assert_eq(bufR->color, color);
 864                 ret = params_->reportHit(
 865                                 rR.fw ? (ebwtFwR?  bufR->patFw  :  bufR->patFwRev) :
 866                                             (ebwtFwR?  bufR->patRc  :  bufR->patRcRev),
 867                                 rR.fw ? (ebwtFwR? &bufR->qual    : &bufR->qualRev) :
 868                                         (ebwtFwR? &bufR->qualRev : &bufR->qual),
 869                                 &bufR->name,
 870                                 bufR->color,
 871                                 colorExEnds,
 872                                 snpPhred,
 873                                 refs_,
 874                                 rmap,
 875                                 ebwtFwR,
 876                                 rR.mms,                       // mismatch positions
 877                                 rR.refcs,                     // reference characters for mms
 878                                 rR.numMms,                    // # mismatches
 879                                 make_pair(first, dnstreamOff),// position
 880                                 make_pair(first, upstreamOff),// mate position
 881                                 rL.fw,                        // mate orientation
 882                                 lenL,                         // mate length
 883                                 make_pair(rR.top, rR.bot),    // arrows
 884                                 tlen,                         // textlen
 885                                 lenR,                         // qlen
 886                                 rR.stratum,                   // alignment stratum
 887                                 rR.cost,                      // cost, including quality penalty
 888                                 oms,                          // # other hits
 889                                 bufR->patid,
 890                                 bufR->seed,
 891                                 pairFw ? 2 : 1);
 892                 return ret;
 893         }
 894
 895         bool report(const Range& rL, // range for upstream mate
 896                     const Range& rR, // range for downstream mate
 897                     uint32_t first,  // ref idx
 898                     uint32_t upstreamOff, // offset for upstream mate
 899                     uint32_t dnstreamOff, // offset for downstream mate
 900                     uint32_t tlen, // length of ref
 901                     bool pairFw,   // whether the pair is being mapped to fw strand
 902                     const ReferenceMap* rmap)
 903         {
 904                 return report(rL, rR, first, upstreamOff,
 905                               dnstreamOff, tlen,
 906                               pairFw, rL.ebwt->fw(), rR.ebwt->fw(), rmap);
 907         }
 908
 909         /**
 910          * Given a vector of reference positions where one of the two mates
 911          * (the "anchor" mate) has aligned, look directly at the reference
 912          * sequence for instances where the other mate (the "outstanding"
 913          * mate) aligns such that mating constraint is satisfied.
 914          *
 915          * This function picks up to 'pick' anchors at random from the
 916          * 'offs' array.  It returns the number that it actually picked.
 917          */
 918         bool resolveOutstandingInRef(const bool off1,
 919                                      const U32Pair& off,
 920                                      const uint32_t tlen,
 921                                      const Range& range)
 922         {
 923                 assert(refs_->loaded());
 924                 assert_lt(off.first, refs_->numRefs());
 925                 // If matchRight is true, then we're trying to align the other
 926                 // mate to the right of the already-aligned mate.  Otherwise,
 927                 // to the left.
 928                 bool matchRight = (off1 ? !doneFw_ : doneFw_);
 929                 // Sequence and quals for mate to be matched
 930                 bool fw = off1 ? fw2_ : fw1_; // whether outstanding mate is fw/rc
 931                 if(doneFw_) fw = !fw;
 932                 // 'seq' gets sequence of outstanding mate w/r/t the forward
 933                 // reference strand
 934                 const String<Dna5>& seq  = fw ? (off1 ? patsrc_->bufb().patFw   :
 935                                                         patsrc_->bufa().patFw)  :
 936                                                 (off1 ? patsrc_->bufb().patRc   :
 937                                                         patsrc_->bufa().patRc);
 938                 // 'seq' gets qualities of outstanding mate w/r/t the forward
 939                 // reference strand
 940                 const String<char>& qual = fw ? (off1 ? patsrc_->bufb().qual  :
 941                                                         patsrc_->bufa().qual) :
 942                                                 (off1 ? patsrc_->bufb().qualRev  :
 943                                                         patsrc_->bufa().qualRev);
 944                 uint32_t qlen = seqan::length(seq);  // length of outstanding mate
 945                 uint32_t alen = (off1 ? patsrc_->bufa().length() :
 946                                         patsrc_->bufb().length());
 947                 int minins = minInsert_;
 948                 int maxins = maxInsert_;
 949                 if(fw1_) {
 950                         minins = max<int>(0, minins - patsrc_->bufa().trimmed5);
 951                         maxins = max<int>(0, maxins - patsrc_->bufa().trimmed5);
 952                 } else {
 953                         minins = max<int>(0, minins - patsrc_->bufa().trimmed3);
 954                         maxins = max<int>(0, maxins - patsrc_->bufa().trimmed3);
 955                 }
 956                 if(fw2_) {
 957                         minins = max<int>(0, minins - patsrc_->bufb().trimmed3);
 958                         maxins = max<int>(0, maxins - patsrc_->bufb().trimmed3);
 959                 } else {
 960                         minins = max<int>(0, minins - patsrc_->bufb().trimmed5);
 961                         maxins = max<int>(0, maxins - patsrc_->bufb().trimmed5);
 962                 }
 963                 assert_geq(minins, 0);
 964                 assert_geq(maxins, 0);
 965                 // Don't even try if either of the mates is longer than the
 966                 // maximum insert size.
 967                 if((uint32_t)maxins <= max(qlen, alen)) {
 968                         return false;
 969                 }
 970                 const uint32_t tidx = off.first;
 971                 const uint32_t toff = off.second;
 972                 // Set begin/end to be a range of all reference
 973                 // positions that are legally permitted to be involved in
 974                 // the alignment of the outstanding mate.  It's up to the
 975                 // callee to worry about how to scan these positions.
 976                 uint32_t begin, end;
 977                 assert_geq(maxins, minins);
 978                 uint32_t insDiff = maxins - minins;
 979                 if(matchRight) {
 980                         end = toff + maxins;
 981                         begin = toff + 1;
 982                         if(qlen < alen) begin += alen-qlen;
 983                         if(end > insDiff + qlen) {
 984                                 begin = max<uint32_t>(begin, end - insDiff - qlen);
 985                         }
 986                         end = min<uint32_t>(refs_->approxLen(tidx), end);
 987                         begin = min<uint32_t>(refs_->approxLen(tidx), begin);
 988                 } else {
 989                         if(toff + alen < (uint32_t)maxins) {
 990                                 begin = 0;
 991                         } else {
 992                                 begin = toff + alen - maxins;
 993                         }
 994                         uint32_t mi = min<uint32_t>(alen, qlen);
 995                         end = toff + mi - 1;
 996                         end = min<uint32_t>(end, toff + alen - minins + qlen - 1);
 997                         if(toff + alen + qlen < (uint32_t)minins + 1) end = 0;
 998                 }
 999                 // Check if there's not enough space in the range to fit an
1000                 // alignment for the outstanding mate.
1001                 if(end - begin < qlen) return false;
1002                 std::vector<Range> ranges;
1003                 std::vector<uint32_t> offs;
1004                 refAligner_->find(1, tidx, refs_, seq, qual, begin, end, ranges,
1005                                   offs, doneFw_ ? &pairs_rc_ : &pairs_fw_,
1006                                   toff, fw);
1007                 assert_eq(ranges.size(), offs.size());
1008                 for(size_t i = 0; i < ranges.size(); i++) {
1009                         Range& r = ranges[i];
1010                         r.fw = fw;
1011                         r.cost |= (r.stratum << 14);
1012                         r.mate1 = !off1;
1013                         const uint32_t result = offs[i];
1014                         // Just copy the known range's top and bot for now
1015                         r.top = range.top;
1016                         r.bot = range.bot;
1017                         bool ebwtLFw = matchRight ? range.ebwt->fw() : true;
1018                         bool ebwtRFw = matchRight ? true : range.ebwt->fw();
1019                         if(report(
1020                                 matchRight ? range : r, // range for upstream mate
1021                                 matchRight ? r : range, // range for downstream mate
1022                                 tidx,                   // ref idx
1023                                 matchRight ? toff : result, // upstream offset
1024                                 matchRight ? result : toff, // downstream offset
1025                                 tlen,       // length of ref
1026                                 !doneFw_,   // whether the pair is being mapped to fw strand
1027                                 ebwtLFw,
1028                                 ebwtRFw,
1029                                 range.ebwt->rmap())) return true;
1030                 }
1031                 return false;
1032         }
1033
1034         /**
1035          * Advance paired-end alignment.
1036          */
1037         void advanceOrientation(bool pairFw) {
1038                 assert(!this->done);
1039                 assert(!*donePair_);
1040                 assert(!*chaseL_ || !*chaseR_);
1041                 if(*chaseL_) {
1042                         assert(!rangeMode_);
1043                         assert(!*delayedchaseL_);
1044                         assert(drL_->foundRange);
1045                         assert(rchase_->foundOff() || rchase_->done);
1046                         if(rchase_->foundOff()) {
1047                                 // Resolve this against the reference loci
1048                                 // determined for the other mate
1049                                 const bool overThresh = (*offsLsz_ + *offsRsz_) > mixedThresh_;
1050                                 if(!this->done && (overThresh || dontReconcile_)) {
1051                                         // Because the total size of both ranges exceeds
1052                                         // our threshold, we're now operating in "mixed
1053                                         // mode"
1054                                         const Range& r = drL_->range();
1055                                         assert(r.repOk());
1056                                         if(verbose_) cout << "Making an attempt to find the outstanding mate" << endl;
1057                                         this->done = resolveOutstandingInRef(
1058                                                         pairFw, rchase_->off(),
1059                                                 r.ebwt->_plen[rchase_->off().first], r);
1060                                         if(++mixedAttempts_ > mixedAttemptLim_) {
1061                                                 // Give up on this pair
1062                                                 *donePair_ = true;
1063                                                 return;
1064                                         }
1065                                 }
1066                                 rchase_->reset();
1067                         } else {
1068                                 assert(rchase_->done);
1069                                 // Forget this range; keep looking for ranges
1070                                 *chaseL_ = false;
1071                                 drL_->foundRange = false;
1072                                 if(verbose_) cout << "Done with chase for first mate" << endl;
1073                                 if(*delayedchaseR_) {
1074                                         // Start chasing the delayed range
1075                                         if(verbose_) cout << "Resuming delayed chase for second mate" << endl;
1076                                         assert(drR_->foundRange);
1077                                         const Range& r = drR_->range();
1078                                         assert(r.repOk());
1079                                         uint32_t top = r.top;
1080                                         uint32_t bot = r.bot;
1081                                         uint32_t qlen = doneFw_? qlen1_ : qlen2_;
1082                                         rchase_->setTopBot(top, bot, qlen, rand_, r.ebwt);
1083                                         *chaseR_ = true;
1084                                         *delayedchaseR_ = false;
1085                                 }
1086                         }
1087                 } else if(*chaseR_) {
1088                         assert(!rangeMode_);
1089                         assert(!*delayedchaseR_);
1090                         assert(drR_->foundRange);
1091                         assert(rchase_->foundOff() || rchase_->done);
1092                         if(rchase_->foundOff()) {
1093                                 // Resolve this against the reference loci
1094                                 // determined for the other mate
1095                                 const bool overThresh = (*offsLsz_ + *offsRsz_) > mixedThresh_;
1096                                 if(!this->done && (overThresh || dontReconcile_)) {
1097                                         // Because the total size of both ranges exceeds
1098                                         // our threshold, we're now operating in "mixed
1099                                         // mode"
1100                                         const Range& r = drR_->range();
1101                                         if(verbose_) cout << "Making an attempt to find the outstanding mate" << endl;
1102                                         this->done = resolveOutstandingInRef(
1103                                                         !pairFw, rchase_->off(),
1104                                                 r.ebwt->_plen[rchase_->off().first], r);
1105                                         if(++mixedAttempts_ > mixedAttemptLim_) {
1106                                                 // Give up on this pair
1107                                                 *donePair_ = true;
1108                                                 return;
1109                                         }
1110                                 }
1111                                 rchase_->reset();
1112                         } else {
1113                                 assert(rchase_->done);
1114                                 // Forget this range; keep looking for ranges
1115                                 *chaseR_ = false;
1116                                 drR_->foundRange = false;
1117                                 if(verbose_) cout << "Done with chase for second mate" << endl;
1118                                 if(*delayedchaseL_) {
1119                                         // Start chasing the delayed range
1120                                         if(verbose_) cout << "Resuming delayed chase for first mate" << endl;
1121                                         assert(drL_->foundRange);
1122                                         const Range& r = drL_->range();
1123                                         assert(r.repOk());
1124                                         uint32_t top = r.top;
1125                                         uint32_t bot = r.bot;
1126                                         uint32_t qlen = doneFw_? qlen2_ : qlen1_;
1127                                         rchase_->setTopBot(top, bot, qlen, rand_, r.ebwt);
1128                                         *chaseL_ = true;
1129                                         *delayedchaseL_ = false;
1130                                 }
1131                         }
1132                 }
1133                 if(!this->done && !*donePair_ && !*chaseL_ && !*chaseR_) {
1134                         // Search for more ranges for whichever mate currently has
1135                         // fewer candidate alignments
1136                         if((*offsLsz_ < *offsRsz_ || drR_->done) && !drL_->done) {
1137                                 // If there are no more ranges for the other mate and
1138                                 // there are no candidate alignments either, then we're
1139                                 // not going to find a paired alignment in this
1140                                 // orientation.
1141                                 if(drR_->done && *offsRsz_ == 0) {
1142                                         // Give up on this orientation
1143                                         if(verbose_) cout << "Giving up on paired orientation " << (pairFw? "fw" : "rc") << " in mate 1" << endl;
1144                                         *donePair_ = true;
1145                                         if(verbose2_) cout << *offsLsz_ << " " << *offsRsz_ << endl;
1146                                         return;
1147                                 }
1148                                 assert(!*delayedchaseL_);
1149                                 if(!drL_->foundRange) drL_->advance(ADV_FOUND_RANGE);
1150                                 if(drL_->foundRange) {
1151 #ifndef NDEBUG
1152                                         {
1153                                                 std::set<int64_t>& s = (pairFw ? allTopsL_fw_ : allTopsL_rc_);
1154                                                 int64_t t = drL_->range().top + 1; // add 1 to avoid 0
1155                                                 if(!drL_->range().ebwt->fw()) t = -t; // invert for bw index
1156                                                 assert(s.find(t) == s.end());
1157                                                 s.insert(t);
1158                                         }
1159 #endif
1160                                         // Add the size of this range to the total for this mate
1161                                         *offsLsz_ += (drL_->range().bot - drL_->range().top);
1162                                         if(*offsRsz_ == 0 && (!dontReconcile_ || *offsLsz_ > 3)) {
1163                                                 // Delay chasing this range; we delay to avoid
1164                                                 // needlessly chasing rows in this range when
1165                                                 // the other mate doesn't end up aligning
1166                                                 // anywhere
1167                                                 if(verbose_) cout << "Delaying a chase for first mate" << endl;
1168                                                 *delayedchaseL_ = true;
1169                                         } else {
1170                                                 // Start chasing this range
1171                                                 if(verbose2_) cout << *offsLsz_ << " " << *offsRsz_ << " " << drL_->range().top << endl;
1172                                                 if(verbose_) cout << "Chasing a range for first mate" << endl;
1173                                                 if(*offsLsz_ > symCeiling_ && *offsRsz_ > symCeiling_) {
1174                                                         // Too many candidates for both mates; abort
1175                                                         // without any more searching
1176                                                         *donePair_ = true;
1177                                                         return;
1178                                                 }
1179                                                 // If this is the first range for both mates,
1180                                                 // choose the smaller range to chase down first
1181                                                 if(*delayedchaseR_ && (*offsRsz_ < *offsLsz_)) {
1182                                                         assert(drR_->foundRange);
1183                                                         *delayedchaseR_ = false;
1184                                                         *delayedchaseL_ = true;
1185                                                         *chaseR_ = true;
1186                                                         const Range& r = drR_->range();
1187                                                         assert(r.repOk());
1188                                                         uint32_t qlen = doneFw_? qlen1_ : qlen2_;
1189                                                         rchase_->setTopBot(r.top, r.bot, qlen, rand_, r.ebwt);
1190                                                 } else {
1191                                                         // Use Burrows-Wheeler for this pair (as
1192                                                         // usual)
1193                                                         *chaseL_ = true;
1194                                                         const Range& r = drL_->range();
1195                                                         uint32_t qlen = doneFw_? qlen2_ : qlen1_;
1196                                                         rchase_->setTopBot(r.top, r.bot, qlen, rand_, r.ebwt);
1197                                                 }
1198                                         }
1199                                 }
1200                         } else if(!drR_->done) {
1201                                 // If there are no more ranges for the other mate and
1202                                 // there are no candidate alignments either, then we're
1203                                 // not going to find a paired alignment in this
1204                                 // orientation.
1205                                 if(drL_->done && *offsLsz_ == 0) {
1206                                         // Give up on this orientation
1207                                         if(verbose_) cout << "Giving up on paired orientation " << (pairFw? "fw" : "rc") << " in mate 2" << endl;
1208                                         if(verbose2_) cout << *offsLsz_ << " " << *offsRsz_ << endl;
1209                                         *donePair_ = true;
1210                                         return;
1211                                 }
1212                                 assert(!*delayedchaseR_);
1213                                 if(!drR_->foundRange) drR_->advance(ADV_FOUND_RANGE);
1214                                 if(drR_->foundRange) {
1215 #ifndef NDEBUG
1216                                         {
1217                                                 std::set<int64_t>& s = (pairFw ? allTopsR_fw_ : allTopsR_rc_);
1218                                                 int64_t t = drR_->range().top + 1; // add 1 to avoid 0
1219                                                 if(!drR_->range().ebwt->fw()) t = -t; // invert for bw index
1220                                                 assert(s.find(t) == s.end());
1221                                                 s.insert(t);
1222                                         }
1223 #endif
1224                                         // Add the size of this range to the total for this mate
1225                                         *offsRsz_ += (drR_->range().bot - drR_->range().top);
1226                                         if(*offsLsz_ == 0 && (!dontReconcile_ || *offsRsz_ > 3)) {
1227                                                 // Delay chasing this range; we delay to avoid
1228                                                 // needlessly chasing rows in this range when
1229                                                 // the other mate doesn't end up aligning
1230                                                 // anywhere
1231                                                 if(verbose_) cout << "Delaying a chase for second mate" << endl;
1232                                                 *delayedchaseR_ = true;
1233                                         } else {
1234                                                 // Start chasing this range
1235                                                 if(verbose2_) cout << *offsLsz_ << " " << *offsRsz_ << " " << drR_->range().top << endl;
1236                                                 if(verbose_) cout << "Chasing a range for second mate" << endl;
1237                                                 if(*offsLsz_ > symCeiling_ && *offsRsz_ > symCeiling_) {
1238                                                         // Too many candidates for both mates; abort
1239                                                         // without any more searching
1240                                                         *donePair_ = true;
1241                                                         return;
1242                                                 }
1243                                                 // If this is the first range for both mates,
1244                                                 // choose the smaller range to chase down first
1245                                                 if(*delayedchaseL_ && *offsLsz_ < *offsRsz_) {
1246                                                         assert(drL_->foundRange);
1247                                                         *delayedchaseL_ = false;
1248                                                         *delayedchaseR_ = true;
1249                                                         *chaseL_ = true;
1250                                                         const Range& r = drL_->range();
1251                                                         assert(r.repOk());
1252                                                         uint32_t qlen = doneFw_? qlen2_ : qlen1_;
1253                                                         rchase_->setTopBot(r.top, r.bot, qlen, rand_, r.ebwt);
1254                                                 } else {
1255                                                         // Use Burrows-Wheeler for this pair (as
1256                                                         // usual)
1257                                                         *chaseR_ = true;
1258                                                         const Range& r = drR_->range();
1259                                                         assert(r.repOk());
1260                                                         uint32_t qlen = doneFw_? qlen1_ : qlen2_;
1261                                                         rchase_->setTopBot(r.top, r.bot, qlen, rand_, r.ebwt);
1262                                                 }
1263                                         }
1264                                 }
1265                         } else {
1266                                 // Finished processing ranges for both mates
1267                                 assert(drL_->done && drR_->done);
1268                                 *donePair_ = true;
1269                         }
1270                 }
1271         }
1272
1273         const BitPairReference* refs_;
1274
1275         PatternSourcePerThread *patsrc_;
1276         uint32_t qlen1_;
1277         uint32_t qlen2_;
1278
1279         // Progress state
1280         bool doneFw_;   // finished with forward orientation of both mates?
1281         bool doneFwFirst_;
1282
1283         bool chase1Fw_;
1284         bool chase1Rc_;
1285         bool chase2Fw_;
1286         bool chase2Rc_;
1287
1288         bool delayedChase1Fw_;
1289         bool delayedChase1Rc_;
1290         bool delayedChase2Fw_;
1291         bool delayedChase2Rc_;
1292
1293         // For searching for outstanding mates
1294         RefAligner<String<Dna5> >* refAligner_;
1295
1296         // Temporary HitSink; to be deleted
1297         const HitSinkPerThreadFactory& sinkPtFactory_;
1298         HitSinkPerThread* sinkPt_;
1299
1300         // State for alignment
1301         EbwtSearchParams<String<Dna> >* params_;
1302
1303         // Paired-end boundaries
1304         const uint32_t minInsert_;
1305         const uint32_t maxInsert_;
1306
1307         // Don't attempt pairwise all-versus-all style of mate
1308         // reconciliation; just rely on mixed mode
1309         const bool dontReconcile_;
1310
1311         // If both mates in a given orientation align >= symCeiling times,
1312         // then immediately give up
1313         const uint32_t symCeiling_;
1314
1315         // If the total number of alignments for both mates in a given
1316         // orientation exceeds mixedThresh, then switch to mixed mode
1317         const uint32_t mixedThresh_;
1318         const uint32_t mixedAttemptLim_;
1319         uint32_t mixedAttempts_;
1320
1321         // Orientation of upstream/downstream mates when aligning to
1322         // forward strand
1323         const bool fw1_;
1324         const bool fw2_;
1325
1326         // State for getting alignments from ranges statefully
1327         RangeChaser<String<Dna> >* rchase_;
1328
1329         // true -> be talkative
1330         bool verbose_;
1331         // true -> suppress warnings
1332         bool quiet_;
1333
1334         int maxBts_;
1335         ChunkPool *pool_;
1336         int *btCnt_;
1337
1338         // Range-finding state for first mate
1339         TDriver*      driver1Fw_;
1340         TDriver*      driver1Rc_;
1341         U32PairVec    offs1FwArr_[32];
1342         TRangeVec     ranges1FwArr_[32];
1343         uint32_t      offs1FwSz_; // total size of all ranges found in this category
1344         U32PairVec    offs1RcArr_[32];
1345         TRangeVec     ranges1RcArr_[32];
1346         uint32_t      offs1RcSz_; // total size of all ranges found in this category
1347
1348         // Range-finding state for second mate
1349         TDriver*      driver2Fw_;
1350         TDriver*      driver2Rc_;
1351         U32PairVec    offs2FwArr_[32];
1352         TRangeVec     ranges2FwArr_[32];
1353         uint32_t      offs2FwSz_; // total size of all ranges found in this category
1354         U32PairVec    offs2RcArr_[32];
1355         TRangeVec     ranges2RcArr_[32];
1356         uint32_t      offs2RcSz_; // total size of all ranges found in this category
1357
1358         bool&       chaseL_fw_;
1359         bool&       chaseR_fw_;
1360         bool&       delayedchaseL_fw_;
1361         bool&       delayedchaseR_fw_;
1362         TDriver&    drL_fw_;
1363         TDriver&    drR_fw_;
1364         U32PairVec* offsLarr_fw_;
1365         U32PairVec* offsRarr_fw_;
1366         TRangeVec*  rangesLarr_fw_;
1367         TRangeVec*  rangesRarr_fw_;
1368         uint32_t&   offsLsz_fw_;
1369         uint32_t&   offsRsz_fw_;
1370
1371         bool&       chaseL_rc_;
1372         bool&       chaseR_rc_;
1373         bool&       delayedchaseL_rc_;
1374         bool&       delayedchaseR_rc_;
1375         TDriver&    drL_rc_;
1376         TDriver&    drR_rc_;
1377         U32PairVec* offsLarr_rc_;
1378         U32PairVec* offsRarr_rc_;
1379         TRangeVec*  rangesLarr_rc_;
1380         TRangeVec*  rangesRarr_rc_;
1381         uint32_t&   offsLsz_rc_;
1382         uint32_t&   offsRsz_rc_;
1383
1384         bool*       chaseL_;
1385         bool*       chaseR_;
1386         bool*       delayedchaseL_;
1387         bool*       delayedchaseR_;
1388         TDriver*    drL_;
1389         TDriver*    drR_;
1390         U32PairVec* offsLarr_;
1391         U32PairVec* offsRarr_;
1392         TRangeVec*  rangesLarr_;
1393         TRangeVec*  rangesRarr_;
1394         uint32_t*   offsLsz_;
1395         uint32_t*   offsRsz_;
1396         bool*       donePair_;
1397         bool        fwL_;
1398         bool        fwR_;
1399
1400         /// For keeping track of paired alignments that have already been
1401         /// found for the forward and reverse-comp pair orientations
1402         TSetPairs   pairs_fw_;
1403         TSetPairs   pairs_rc_;
1404
1405 #ifndef NDEBUG
1406         std::set<int64_t> allTopsL_fw_;
1407         std::set<int64_t> allTopsR_fw_;
1408         std::set<int64_t> allTopsL_rc_;
1409         std::set<int64_t> allTopsR_rc_;
1410 #endif
1411
1412         bool verbose2_;
1413 };
1414
1415 /**
1416  * Helper struct that holds a Range together with the coordinates where it al
1417  */
1418 struct RangeWithCoords {
1419         Range r;
1420         U32Pair h;
1421 };
1422
1423 /**
1424  * An aligner for finding paired alignments while operating entirely
1425  * within the Burrows-Wheeler domain.
1426  */
1427 template<typename TRangeSource>
1428 class PairedBWAlignerV2 : public Aligner {
1429
1430         typedef std::pair<uint32_t,uint32_t> U32Pair;
1431         typedef std::vector<U32Pair> U32PairVec;
1432         typedef std::vector<Range> TRangeVec;
1433         typedef RangeSourceDriver<TRangeSource> TDriver;
1434         typedef std::pair<uint64_t, uint64_t> TU64Pair;
1435         typedef std::set<TU64Pair> TSetPairs;
1436
1437 public:
1438         PairedBWAlignerV2(
1439                 EbwtSearchParams<String<Dna> >* params,
1440                 EbwtSearchParams<String<Dna> >* paramsSe1,
1441                 EbwtSearchParams<String<Dna> >* paramsSe2,
1442                 TDriver* driver,
1443                 RefAligner<String<Dna5> >* refAligner,
1444                 RangeChaser<String<Dna> >* rchase,
1445                 HitSink& sink,
1446                 const HitSinkPerThreadFactory& sinkPtFactory,
1447                 HitSinkPerThread* sinkPt,
1448                 HitSinkPerThread* sinkPtSe1,
1449                 HitSinkPerThread* sinkPtSe2,
1450                 bool fw1, bool fw2,
1451                 uint32_t minInsert,
1452                 uint32_t maxInsert,
1453                 uint32_t mixedAttemptLim,
1454                 const BitPairReference* refs,
1455                 bool rangeMode,
1456                 bool verbose,
1457                 bool quiet,
1458                 int maxBts,
1459                 ChunkPool *pool,
1460                 int *btCnt) :
1461                 Aligner(true, rangeMode),
1462                 refs_(refs),
1463                 patsrc_(NULL),
1464                 qlen1_(0), qlen2_(0),
1465                 chase_(false),
1466                 donePe_(false),
1467                 doneSe1_(false),
1468                 doneSe2_(false),
1469                 refAligner_(refAligner),
1470                 sinkPtFactory_(sinkPtFactory),
1471                 sinkPt_(sinkPt),
1472                 sinkPtSe1_(sinkPtSe1),
1473                 sinkPtSe2_(sinkPtSe2),
1474                 params_(params),
1475                 paramsSe1_(paramsSe1),
1476                 paramsSe2_(paramsSe2),
1477                 minInsert_(minInsert),
1478                 maxInsert_(maxInsert),
1479                 mixedAttemptLim_(mixedAttemptLim),
1480                 mixedAttempts_(0),
1481                 fw1_(fw1), fw2_(fw2),
1482                 rchase_(rchase),
1483                 driver_(driver),
1484                 pool_(pool),
1485                 verbose_(verbose),
1486                 quiet_(quiet),
1487                 maxBts_(maxBts),
1488                 btCnt_(btCnt)
1489         {
1490                 assert(sinkPt_ != NULL);
1491                 assert(params_ != NULL);
1492                 assert(driver_ != NULL);
1493         }
1494
1495         virtual ~PairedBWAlignerV2() {
1496                 delete driver_; driver_ = NULL;
1497                 delete params_; params_ = NULL;
1498                 if(paramsSe1_ != NULL) {
1499                         delete paramsSe1_; paramsSe1_ = NULL;
1500                         delete paramsSe2_; paramsSe2_ = NULL;
1501                 }
1502                 delete rchase_; rchase_ = NULL;
1503                 delete[] btCnt_; btCnt_ = NULL;
1504                 delete refAligner_; refAligner_ = NULL;
1505                 sinkPtFactory_.destroy(sinkPt_); sinkPt_ = NULL;
1506                 if(sinkPtSe1_ != NULL) {
1507                         sinkPtFactory_.destroy(sinkPtSe1_); sinkPtSe1_ = NULL;
1508                         sinkPtFactory_.destroy(sinkPtSe2_); sinkPtSe2_ = NULL;
1509                 }
1510         }
1511
1512         /**
1513          * Prepare this aligner for the next read.
1514          */
1515         virtual void setQuery(PatternSourcePerThread* patsrc) {
1516                 assert(!patsrc->bufa().empty());
1517                 Aligner::setQuery(patsrc); // set fields & random seed
1518                 assert(!patsrc->bufb().empty());
1519                 // Give all of the drivers pointers to the relevant read info
1520                 patsrc_ = patsrc;
1521                 pool_->reset(&patsrc->bufa().name, patsrc->patid());
1522                 if(patsrc->bufa().length() < 4 || patsrc->bufb().length() < 4) {
1523                         if(!quiet_) {
1524                                 cerr << "Warning: Skipping pair " << patsrc->bufa().name
1525                                      << " because a mate is less than 4 characters long" << endl;
1526                         }
1527                         this->done = true;
1528                         sinkPt_->finishRead(*patsrc_, true, true);
1529                         return;
1530                 }
1531                 driver_->setQuery(patsrc, NULL);
1532                 qlen1_ = patsrc_->bufa().length();
1533                 qlen2_ = patsrc_->bufb().length();
1534                 if(btCnt_ != NULL) (*btCnt_) = maxBts_;
1535                 mixedAttempts_ = 0;
1536                 // Neither orientation is done
1537                 this->done = false;
1538                 // No ranges are being chased yet
1539                 chase_ = false;
1540                 donePe_ = doneSe1_ = doneSe2_ = false;
1541                 pairs_fw_.clear();
1542                 pairs_rc_.clear();
1543         }
1544
1545         /**
1546          * Advance the aligner by one memory op.  Return true iff we're
1547          * done with this read.
1548          *
1549          * A call to this function does one of many things:
1550          * 1. Advance a RangeSourceDriver and check if it found a new range
1551          * 2. Advance a RowChaseDriver and check if it found a reference
1552          *    offset for a an alignment in a range
1553          */
1554         virtual bool advance() {
1555                 assert(!this->done);
1556                 if(chase_) {
1557                         assert(!rangeMode_); // chasing ranges
1558                         if(!rchase_->foundOff() && !rchase_->done) {
1559                                 rchase_->advance();
1560                                 return false;
1561                         }
1562                         assert(rchase_->foundOff() || rchase_->done);
1563                         if(rchase_->foundOff()) {
1564                                 const Range& r = driver_->range();
1565                                 assert(r.repOk());
1566                                 resolveOutstanding(
1567                                         rchase_->off(),
1568                                         r.ebwt->_plen[rchase_->off().first], r);
1569                                 rchase_->reset();
1570                         } else {
1571                                 assert(rchase_->done);
1572                                 // Forget this range; keep looking for ranges
1573                                 chase_ = false;
1574                                 this->done = driver_->done;
1575                         }
1576                 }
1577
1578                 if(!this->done && !chase_) {
1579                         // Search for more ranges for whichever mate currently has
1580                         // fewer candidate alignments
1581                         if(!driver_->done) {
1582                                 if(!this->done) {
1583                                         //
1584                                         // Check whether any of the PE/SE possibilities
1585                                         // have become impossible due to the minCost
1586                                         //
1587                                         if(!donePe_) {
1588                                                 assert(!this->done);
1589                                                 donePe_ = sinkPt_->irrelevantCost(driver_->minCost);
1590                                                 if(donePe_ && (!sinkPt_->empty() || sinkPtSe1_ == NULL)) {
1591                                                         // Paired-end alignment(s) were found, no
1592                                                         // more will be found, and no unpaired
1593                                                         // alignments are requested, so stop
1594                                                         this->done = true;
1595                                                 }
1596                                                 if(donePe_ && sinkPtSe1_ != NULL) {
1597                                                         // Note: removeMate affects minCost
1598                                                         if(doneSe1_) driver_->removeMate(1);
1599                                                         if(doneSe2_) driver_->removeMate(2);
1600                                                 }
1601                                         }
1602                                         if(!this->done && sinkPtSe1_ != NULL) {
1603                                                 if(!doneSe1_) {
1604                                                         doneSe1_ = sinkPtSe1_->irrelevantCost(driver_->minCost);
1605                                                         if(doneSe1_ && donePe_) driver_->removeMate(1);
1606                                                 }
1607                                                 if(!doneSe2_) {
1608                                                         doneSe2_ = sinkPtSe2_->irrelevantCost(driver_->minCost);
1609                                                         if(doneSe2_ && donePe_) driver_->removeMate(2);
1610                                                 }
1611                                                 // Do Se1 again, because removing Se2 may have
1612                                                 // nudged minCost over the threshold
1613                                                 if(!doneSe1_) {
1614                                                         doneSe1_ = sinkPtSe1_->irrelevantCost(driver_->minCost);
1615                                                         if(doneSe1_ && donePe_) driver_->removeMate(1);
1616                                                 }
1617                                                 if(doneSe1_ && doneSe2_) assert(donePe_);
1618                                                 this->done = donePe_ && doneSe1_ && doneSe2_;
1619                                         }
1620
1621                                         if(!this->done) {
1622                                                 if(sinkPtSe1_ != NULL) {
1623                                                         assert(doneSe1_ || !sinkPtSe1_->irrelevantCost(driver_->minCost));
1624                                                         assert(doneSe2_ || !sinkPtSe2_->irrelevantCost(driver_->minCost));
1625                                                 }
1626                                                 assert(donePe_ || !sinkPt_->irrelevantCost(driver_->minCost));
1627                                                 driver_->advance(ADV_COST_CHANGES);
1628                                         }
1629                                 }
1630                                 if(driver_->foundRange) {
1631                                         // Use Burrows-Wheeler for this pair (as usual)
1632                                         chase_ = true;
1633                                         driver_->foundRange = false;
1634                                         const Range& r = driver_->range();
1635                                         assert(r.repOk());
1636                                         rchase_->setTopBot(r.top, r.bot,
1637                                                            r.mate1 ? qlen1_ : qlen2_,
1638                                                            rand_, r.ebwt);
1639                                 }
1640                         } else {
1641                                 this->done = true;
1642                         }
1643                 }
1644
1645                 if(this->done) {
1646                         bool reportedPe = (sinkPt_->finishRead(*patsrc_, true, true) > 0);
1647                         if(sinkPtSe1_ != NULL) {
1648                                 sinkPtSe1_->finishRead(*patsrc_, !reportedPe, false);
1649                                 sinkPtSe2_->finishRead(*patsrc_, !reportedPe, false);
1650                         }
1651                 }
1652                 return this->done;
1653         }
1654
1655 protected:
1656
1657         /**
1658          * Helper for reporting a pair of alignments.  As of now, we report
1659          * a paired alignment by reporting two consecutive alignments, one
1660          * for each mate.
1661          */
1662         bool report(const Range& rL, // range for upstream mate
1663                     const Range& rR, // range for downstream mate
1664                     uint32_t first,  // ref idx
1665                     uint32_t upstreamOff, // offset for upstream mate
1666                     uint32_t dnstreamOff, // offset for downstream mate
1667                     uint32_t tlen, // length of ref
1668                     bool pairFw,   // whether the pair is being mapped to fw strand
1669                     bool ebwtFwL,
1670                     bool ebwtFwR,
1671                     const ReferenceMap *rmap)
1672         {
1673                 assert_lt(upstreamOff, dnstreamOff);
1674                 uint32_t spreadL = rL.bot - rL.top;
1675                 uint32_t spreadR = rR.bot - rR.top;
1676                 uint32_t oms = min(spreadL, spreadR) - 1;
1677                 ReadBuf* bufL = pairFw ? bufa_ : bufb_;
1678                 ReadBuf* bufR = pairFw ? bufb_ : bufa_;
1679                 uint32_t lenL = pairFw ? alen_ : blen_;
1680                 uint32_t lenR = pairFw ? blen_ : alen_;
1681                 bool ret;
1682                 assert(!params_->sink().exceededOverThresh());
1683                 params_->setFw(rL.fw);
1684                 assert_eq(bufL->color, color);
1685                 // Print upstream mate first
1686                 ret = params_->reportHit(
1687                                 rL.fw ? (ebwtFwL?  bufL->patFw  :  bufL->patFwRev) :
1688                                         (ebwtFwL?  bufL->patRc  :  bufL->patRcRev),
1689                                 rL.fw ? (ebwtFwL? &bufL->qual    : &bufL->qualRev) :
1690                                         (ebwtFwL? &bufL->qualRev : &bufL->qual),
1691                                 &bufL->name,
1692                                 bufL->color,
1693                                 colorExEnds,
1694                                 snpPhred,
1695                                 refs_,
1696                                 rmap,
1697                                 ebwtFwL,
1698                                 rL.mms,                       // mismatch positions
1699                                 rL.refcs,                     // reference characters for mms
1700                                 rL.numMms,                    // # mismatches
1701                                 make_pair(first, upstreamOff),// position
1702                                 make_pair(first, dnstreamOff),// mate position
1703                                 rR.fw,                        // mate orientation
1704                                 lenR,                         // mate length
1705                                 make_pair(rL.top, rL.bot),    // arrows
1706                                 tlen,                         // textlen
1707                                 lenL,                         // qlen
1708                                 rL.stratum,                   // alignment stratum
1709                                 rL.cost,                      // cost, including quality penalty
1710                                 oms,                          // # other hits
1711                                 bufL->patid,
1712                                 bufL->seed,
1713                                 pairFw ? 1 : 2);
1714                 if(ret) {
1715                         return true; // can happen when -m is set
1716                 }
1717                 params_->setFw(rR.fw);
1718                 assert_eq(bufR->color, color);
1719                 ret = params_->reportHit(
1720                                 rR.fw ? (ebwtFwR?  bufR->patFw  :  bufR->patFwRev) :
1721                                         (ebwtFwR?  bufR->patRc  :  bufR->patRcRev),
1722                                 rR.fw ? (ebwtFwR? &bufR->qual    : &bufR->qualRev) :
1723                                         (ebwtFwR? &bufR->qualRev : &bufR->qual),
1724                                 &bufR->name,
1725                                 bufR->color,
1726                                 colorExEnds,
1727                                 snpPhred,
1728                                 refs_,
1729                                 rmap,
1730                                 ebwtFwR,
1731                                 rR.mms,                       // mismatch positions
1732                                 rR.refcs,                     // reference characters for mms
1733                                 rR.numMms,                    // # mismatches
1734                                 make_pair(first, dnstreamOff),// position
1735                                 make_pair(first, upstreamOff),// mate position
1736                                 rL.fw,                        // mate orientation
1737                                 lenL,                         // mate length
1738                                 make_pair(rR.top, rR.bot),    // arrows
1739                                 tlen,                         // textlen
1740                                 lenR,                         // qlen
1741                                 rR.stratum,                   // alignment stratum
1742                                 rR.cost,                      // cost, including quality penalty
1743                                 oms,                          // # other hits
1744                                 bufR->patid,
1745                                 bufR->seed,
1746                                 pairFw ? 2 : 1);
1747                 return ret;
1748         }
1749
1750         /**
1751          * Helper for reporting a pair of alignments.  As of now, we report
1752          * a paired alignment by reporting two consecutive alignments, one
1753          * for each mate.
1754          */
1755         void reportSe(const Range& r, U32Pair h, uint32_t tlen) {
1756                 EbwtSearchParams<String<Dna> >*params = (r.mate1 ? paramsSe1_ : paramsSe2_);
1757                 assert(!(r.mate1 ? doneSe1_ : doneSe2_));
1758                 params->setFw(r.fw);
1759                 ReadBuf* buf = r.mate1 ? bufa_ : bufb_;
1760                 bool ebwtFw = r.ebwt->fw();
1761                 uint32_t len = r.mate1 ? alen_ : blen_;
1762                 assert_eq(buf->color, color);
1763                 // Print upstream mate first
1764                 if(params->reportHit(
1765                         r.fw ? (ebwtFw?  buf->patFw   :  buf->patFwRev) :
1766                                (ebwtFw?  buf->patRc   :  buf->patRcRev),
1767                         r.fw ? (ebwtFw? &buf->qual    : &buf->qualRev) :
1768                                (ebwtFw? &buf->qualRev : &buf->qual),
1769                         &buf->name,
1770                         buf->color,
1771                         colorExEnds,
1772                         snpPhred,
1773                         refs_,
1774                         r.ebwt->rmap(),
1775                         ebwtFw,
1776                         r.mms,                   // mismatch positions
1777                         r.refcs,                 // reference characters for mms
1778                         r.numMms,                // # mismatches
1779                         h,                       // position
1780                         make_pair(0, 0),         // (bogus) mate coords
1781                         true,                    // (bogus) mate orientation
1782                         0,                       // (bogus) mate length
1783                         make_pair(r.top, r.bot), // arrows
1784                         tlen,                    // textlen
1785                         len,                     // qlen
1786                         r.stratum,               // alignment stratum
1787                         r.cost,                  // cost, including quality penalty
1788                         r.bot - r.top - 1,       // # other hits
1789                         buf->patid,
1790                         buf->seed,
1791                         0))
1792                 {
1793                         if(r.mate1) doneSe1_ = true;
1794                         else        doneSe2_ = true;
1795                         if(donePe_) driver_->removeMate(r.mate1 ? 1 : 2);
1796                 }
1797         }
1798
1799         void resolveOutstanding(const U32Pair& off,
1800                                 const uint32_t tlen,
1801                                 const Range& range)
1802         {
1803                 assert(!this->done);
1804                 if(!donePe_) {
1805                         bool ret = resolveOutstandingInRef(off, tlen, range);
1806                         if(++mixedAttempts_ > mixedAttemptLim_ || ret) {
1807                                 // Give up on this pair
1808                                 donePe_ = true;
1809                                 if(sinkPtSe1_ != NULL) {
1810                                         if(doneSe1_) driver_->removeMate(1);
1811                                         if(doneSe2_) driver_->removeMate(2);
1812                                 }
1813                         }
1814                         this->done = (donePe_ && (!sinkPt_->empty() || sinkPtSe1_ == NULL || (doneSe1_ && doneSe2_)));
1815                 }
1816                 if(!this->done && sinkPtSe1_ != NULL) {
1817                         bool doneSe = (range.mate1 ? doneSe1_ : doneSe2_);
1818                         if(!doneSe) {
1819                                 // Hold onto this single-end alignment in case we don't
1820                                 // find any paired alignments
1821                                 reportSe(range, off, tlen);
1822                         }
1823                         this->done = doneSe1_ && doneSe2_ && donePe_;
1824                 }
1825         }
1826
1827         /**
1828          * Given a vector of reference positions where one of the two mates
1829          * (the "anchor" mate) has aligned, look directly at the reference
1830          * sequence for instances where the other mate (the "outstanding"
1831          * mate) aligns such that mating constraint is satisfied.
1832          *
1833          * This function picks up to 'pick' anchors at random from the
1834          * 'offs' array.  It returns the number that it actually picked.
1835          */
1836         bool resolveOutstandingInRef(const U32Pair& off,
1837                                      const uint32_t tlen,
1838                                      const Range& range)
1839         {
1840                 assert(!donePe_);
1841                 assert(refs_->loaded());
1842                 assert_lt(off.first, refs_->numRefs());
1843                 // pairFw = true if the anchor indicates that the pair will
1844                 // align in its forward orientation (i.e. with mate1 to the
1845                 // left of mate2)
1846                 bool pairFw = (range.mate1)? (range.fw == fw1_) : (range.fw == fw2_);
1847                 // matchRight = true, if the opposite mate will be to the right
1848                 // of the anchor mate
1849                 bool matchRight = (pairFw ? range.mate1 : !range.mate1);
1850                 // fw = orientation of the opposite mate
1851                 bool fw = range.mate1 ? fw2_ : fw1_; // whether outstanding mate is fw/rc
1852                 if(!pairFw) fw = !fw;
1853                 // 'seq' = sequence for opposite mate
1854                 const String<Dna5>& seq  =
1855                         fw ? (range.mate1 ? patsrc_->bufb().patFw   :
1856                                         patsrc_->bufa().patFw)  :
1857                          (range.mate1 ? patsrc_->bufb().patRc   :
1858                                         patsrc_->bufa().patRc);
1859                 // 'qual' = qualities for opposite mate
1860                 const String<char>& qual =
1861                         fw ? (range.mate1 ? patsrc_->bufb().qual  :
1862                                             patsrc_->bufa().qual) :
1863                              (range.mate1 ? patsrc_->bufb().qualRev :
1864                                             patsrc_->bufa().qualRev);
1865                 uint32_t qlen = seqan::length(seq);  // length of outstanding mate
1866                 uint32_t alen = (range.mate1 ? patsrc_->bufa().length() :
1867                                                patsrc_->bufb().length());
1868                 int minins = minInsert_;
1869                 int maxins = maxInsert_;
1870                 if(fw1_) {
1871                         minins = max<int>(0, minins - patsrc_->bufa().trimmed5);
1872                         maxins = max<int>(0, maxins - patsrc_->bufa().trimmed5);
1873                 } else {
1874                         minins = max<int>(0, minins - patsrc_->bufa().trimmed3);
1875                         maxins = max<int>(0, maxins - patsrc_->bufa().trimmed3);
1876                 }
1877                 if(fw2_) {
1878                         minins = max<int>(0, minins - patsrc_->bufb().trimmed3);
1879                         maxins = max<int>(0, maxins - patsrc_->bufb().trimmed3);
1880                 } else {
1881                         minins = max<int>(0, minins - patsrc_->bufb().trimmed5);
1882                         maxins = max<int>(0, maxins - patsrc_->bufb().trimmed5);
1883                 }
1884                 assert_geq(minins, 0);
1885                 assert_geq(maxins, 0);
1886                 // Don't even try if either of the mates is longer than the
1887                 // maximum insert size.
1888                 if((uint32_t)maxins <= max(qlen, alen)) {
1889                         return false;
1890                 }
1891                 const uint32_t tidx = off.first;  // text id where anchor mate hit
1892                 const uint32_t toff = off.second; // offset where anchor mate hit
1893                 // Set begin/end to the range of reference positions where
1894                 // outstanding mate may align while fulfilling insert-length
1895                 // constraints.
1896                 uint32_t begin, end;
1897                 assert_geq(maxins, minins);
1898                 uint32_t insDiff = maxins - minins;
1899                 if(matchRight) {
1900                         end = toff + maxins;
1901                         begin = toff + 1;
1902                         if(qlen < alen) begin += alen-qlen;
1903                         if(end > insDiff + qlen) {
1904                                 begin = max<uint32_t>(begin, end - insDiff - qlen);
1905                         }
1906                         end = min<uint32_t>(refs_->approxLen(tidx), end);
1907                         begin = min<uint32_t>(refs_->approxLen(tidx), begin);
1908                 } else {
1909                         if(toff + alen < (uint32_t)maxins) {
1910                                 begin = 0;
1911                         } else {
1912                                 begin = toff + alen - maxins;
1913                         }
1914                         uint32_t mi = min<uint32_t>(alen, qlen);
1915                         end = toff + mi - 1;
1916                         end = min<uint32_t>(end, toff + alen - minins + qlen - 1);
1917                         if(toff + alen + qlen < (uint32_t)(minins + 1)) end = 0;
1918                 }
1919                 // Check if there's not enough space in the range to fit an
1920                 // alignment for the outstanding mate.
1921                 if(end - begin < qlen) return false;
1922                 std::vector<Range> ranges;
1923                 std::vector<uint32_t> offs;
1924                 refAligner_->find(1, tidx, refs_, seq, qual, begin, end, ranges,
1925                                   offs, pairFw ? &pairs_fw_ : &pairs_rc_,
1926                                   toff, fw);
1927                 assert_eq(ranges.size(), offs.size());
1928                 for(size_t i = 0; i < ranges.size(); i++) {
1929                         Range& r = ranges[i];
1930                         r.fw = fw;
1931                         r.cost |= (r.stratum << 14);
1932                         r.mate1 = !range.mate1;
1933                         const uint32_t result = offs[i];
1934                         // Just copy the known range's top and bot for now
1935                         r.top = range.top;
1936                         r.bot = range.bot;
1937                         bool ebwtLFw = matchRight ? range.ebwt->fw() : true;
1938                         bool ebwtRFw = matchRight ? true : range.ebwt->fw();
1939                         if(report(
1940                                 matchRight ? range : r, // range for upstream mate
1941                                 matchRight ? r : range, // range for downstream mate
1942                                 tidx,                   // ref idx
1943                                 matchRight ? toff : result, // upstream offset
1944                                 matchRight ? result : toff, // downstream offset
1945                                 tlen,       // length of ref
1946                                 pairFw,     // whether the pair is being mapped to fw strand
1947                                 ebwtLFw,
1948                                 ebwtRFw,
1949                                 range.ebwt->rmap())) return true;
1950                 }
1951                 return false;
1952         }
1953
1954         const BitPairReference* refs_;
1955
1956         PatternSourcePerThread *patsrc_;
1957         uint32_t qlen1_, qlen2_;
1958         bool chase_;
1959
1960         // true -> we're no longer shooting for paired-end alignments;
1961         // just collecting single-end ones
1962         bool donePe_, doneSe1_, doneSe2_;
1963
1964         // For searching for outstanding mates
1965         RefAligner<String<Dna5> >* refAligner_;
1966
1967         // Temporary HitSink; to be deleted
1968         const HitSinkPerThreadFactory& sinkPtFactory_;
1969         HitSinkPerThread* sinkPt_;
1970         HitSinkPerThread* sinkPtSe1_, * sinkPtSe2_;
1971
1972         // State for alignment
1973         EbwtSearchParams<String<Dna> >* params_;
1974         // for single-end:
1975         EbwtSearchParams<String<Dna> >* paramsSe1_, * paramsSe2_;
1976
1977         // Paired-end boundaries
1978         const uint32_t minInsert_;
1979         const uint32_t maxInsert_;
1980
1981         const uint32_t mixedAttemptLim_;
1982         uint32_t mixedAttempts_;
1983
1984         // Orientation of upstream/downstream mates when aligning to
1985         // forward strand
1986         const bool fw1_, fw2_;
1987
1988         // State for getting alignments from ranges statefully
1989         RangeChaser<String<Dna> >* rchase_;
1990
1991         // Range-finding state for first mate
1992         TDriver* driver_;
1993
1994         // Pool for distributing chunks of best-first path descriptor memory
1995         ChunkPool *pool_;
1996
1997         bool verbose_;
1998         bool quiet_;
1999
2000         int maxBts_; // maximum allowed # backtracks
2001         int *btCnt_; // current backtrack count
2002
2003         /// For keeping track of paired alignments that have already been
2004         /// found for the forward and reverse-comp pair orientations
2005         TSetPairs pairs_fw_, pairs_rc_;
2006 };
2007
2008 #endif /* ALIGNER_H_ */