src/toolkit/bamtools_resolve.cpp

   1 // ***************************************************************************
   2 // bamtools_resolve.cpp (c) 2011
   3 // Marth Lab, Department of Biology, Boston College
   4 // All rights reserved.
   5 // ---------------------------------------------------------------------------
   6 // Last modified: 28 June 2011
   7 // ---------------------------------------------------------------------------
   8 // Resolves paired-end reads (marking the IsProperPair flag as needed).
   9 // ***************************************************************************
  10
  11 #include "bamtools_resolve.h"
  12 #include "bamtools_version.h"
  13 #include <api/BamReader.h>
  14 #include <api/BamWriter.h>
  15 #include <utils/bamtools_options.h>
  16 #include <utils/bamtools_utilities.h>
  17 using namespace BamTools;
  18
  19 #include <algorithm>
  20 #include <cassert>
  21 #include <cctype>
  22 #include <cstdio>
  23 #include <cstdlib>
  24 #include <fstream>
  25 #include <iostream>
  26 #include <map>
  27 #include <sstream>
  28 #include <string>
  29 #include <utility>
  30 #include <vector>
  31 using namespace std;
  32
  33 // --------------------------------------------------------------------------
  34 // general ResolveTool constants
  35 // --------------------------------------------------------------------------
  36
  37 static const int NUM_MODELS = 8;
  38 static const string READ_GROUP_TAG = "RG";
  39 static const double DEFAULT_CONFIDENCE_INTERVAL = 0.9973;
  40 static const double DEFAULT_UNUSEDMODEL_THRESHOLD = 0.1;
  41
  42 // --------------------------------------------------------------------------
  43 // stats file constants
  44 // --------------------------------------------------------------------------
  45
  46 // basic char/string constants
  47 static const char COMMENT_CHAR     = '#';
  48 static const char OPEN_BRACE_CHAR  = '[';
  49 static const char CLOSE_BRACE_CHAR = ']';
  50 static const char EQUAL_CHAR       = '=';
  51 static const char TAB_CHAR         = '\t';
  52
  53 static const string WHITESPACE_CHARS = " \t\n";
  54 static const string TRUE_KEYWORD     = "true";
  55 static const string FALSE_KEYWORD    = "false";
  56
  57 // field counts
  58 static const size_t NUM_OPTIONS_FIELDS    = 2;
  59 static const size_t NUM_READGROUPS_FIELDS = 7;
  60
  61 // header strings
  62 static const string INPUT_TOKEN      = "[Input]";
  63 static const string OPTIONS_TOKEN    = "[Options]";
  64 static const string READGROUPS_TOKEN = "[ReadGroups]";
  65
  66 // option keywords
  67 static const string OPTION_CONFIDENCEINTERVAL   = "ConfidenceInterval";
  68 static const string OPTION_UNUSEDMODELTHRESHOLD = "UnusedModelThreshold";
  69 static const string OPTION_FORCEMARKREADGROUPS  = "ForceMarkReadGroups";
  70
  71 // other string constants
  72 static const string RG_FIELD_DESCRIPTION =
  73     "#<name> <medianFL> <minFL> <maxFL> <topModelID> <nextTopModelID> <isAmbiguous?>";
  74
  75 // --------------------------------------------------------------------------
  76 // unique readname file constants
  77 // --------------------------------------------------------------------------
  78
  79 static const string READNAME_FILE_SUFFIX = ".uniq_names.txt";
  80 static const string DEFAULT_READNAME_FILE = "bt_resolve_TEMP" + READNAME_FILE_SUFFIX;
  81
  82 // --------------------------------------------------------------------------
  83 // ModelType implementation
  84
  85 struct ModelType {
  86
  87     // data members
  88     uint16_t ID;
  89     vector<int32_t> FragmentLengths;
  90
  91     // ctor
  92     ModelType(const uint16_t id)
  93         : ID(id)
  94     {
  95         // preallocate space for 10K fragments per model type
  96         FragmentLengths.reserve(10000);
  97     }
  98
  99     // convenience access to internal fragment lengths vector
 100     vector<int32_t>::iterator begin(void) { return FragmentLengths.begin(); }
 101     vector<int32_t>::const_iterator begin(void) const { return FragmentLengths.begin(); }
 102     void clear(void) { FragmentLengths.clear(); }
 103     vector<int32_t>::iterator end(void) { return FragmentLengths.end(); }
 104     vector<int32_t>::const_iterator end(void) const { return FragmentLengths.end(); }
 105     void push_back(const int32_t& x) { FragmentLengths.push_back(x); }
 106     size_t size(void) const { return FragmentLengths.size(); }
 107
 108     // constants
 109     static const uint16_t DUMMY_ID;
 110 };
 111
 112 const uint16_t ModelType::DUMMY_ID = 100;
 113
 114 bool operator>(const ModelType& lhs, const ModelType& rhs) {
 115     return lhs.size() > rhs.size();
 116 }
 117
 118 uint16_t CalculateModelType(const BamAlignment& al) {
 119
 120     // localize alignment's mate positions & orientations for convenience
 121     const int32_t m1_begin = ( al.IsFirstMate() ? al.Position : al.MatePosition );
 122     const int32_t m2_begin = ( al.IsFirstMate() ? al.MatePosition : al.Position );
 123     const bool m1_isReverseStrand = ( al.IsFirstMate() ? al.IsReverseStrand() : al.IsMateReverseStrand() );
 124     const bool m2_isReverseStrand = ( al.IsFirstMate() ? al.IsMateReverseStrand() : al.IsReverseStrand() );
 125
 126     // determine 'model type'
 127     if ( m1_begin < m2_begin ) {
 128         if ( !m1_isReverseStrand && !m2_isReverseStrand ) return 0; // ID: 1
 129         if ( !m1_isReverseStrand &&  m2_isReverseStrand ) return 1; // ID: 2
 130         if (  m1_isReverseStrand && !m2_isReverseStrand ) return 2; // ID: 3
 131         if (  m1_isReverseStrand &&  m2_isReverseStrand ) return 3; // ID: 4
 132     } else {
 133         if ( !m2_isReverseStrand && !m1_isReverseStrand ) return 4; // ID: 5
 134         if ( !m2_isReverseStrand &&  m1_isReverseStrand ) return 5; // ID: 6
 135         if (  m2_isReverseStrand && !m1_isReverseStrand ) return 6; // ID: 7
 136         if (  m2_isReverseStrand &&  m1_isReverseStrand ) return 7; // ID: 8
 137     }
 138
 139     // unknown model
 140     return ModelType::DUMMY_ID;
 141 }
 142
 143 // --------------------------------------------------------------------------
 144 // ReadGroupResolver implementation
 145
 146 struct ReadGroupResolver {
 147
 148     // data members
 149     int32_t MinFragmentLength;
 150     int32_t MedianFragmentLength;
 151     int32_t MaxFragmentLength;
 152     uint16_t TopModelId;
 153     uint16_t NextTopModelId;
 154     bool IsAmbiguous;
 155     bool HasData;
 156     vector<ModelType> Models;
 157     map<string, bool> ReadNames;
 158
 159     // ctor
 160     ReadGroupResolver(void);
 161
 162     // resolving methods
 163     bool IsValidInsertSize(const BamAlignment& al) const;
 164     bool IsValidOrientation(const BamAlignment& al) const;
 165
 166     // select 2 best models based on observed data
 167     void DetermineTopModels(const string& readGroupName);
 168
 169     // static settings
 170     static double ConfidenceInterval;
 171     static double UnusedModelThreshold;
 172     static void SetConfidenceInterval(const double& ci);
 173     static void SetUnusedModelThreshold(const double& umt);
 174 };
 175
 176 double ReadGroupResolver::ConfidenceInterval   = DEFAULT_CONFIDENCE_INTERVAL;
 177 double ReadGroupResolver::UnusedModelThreshold = DEFAULT_UNUSEDMODEL_THRESHOLD;
 178
 179 ReadGroupResolver::ReadGroupResolver(void)
 180     : MinFragmentLength(0)
 181     , MedianFragmentLength(0)
 182     , MaxFragmentLength(0)
 183     , TopModelId(ModelType::DUMMY_ID)
 184     , NextTopModelId(ModelType::DUMMY_ID)
 185     , IsAmbiguous(false)
 186     , HasData(false)
 187 {
 188     // pre-allocate space for 8 models
 189     Models.reserve(NUM_MODELS);
 190     for ( uint16_t i = 0; i < NUM_MODELS; ++i )
 191         Models.push_back( ModelType(i+1) );
 192 }
 193
 194 bool ReadGroupResolver::IsValidInsertSize(const BamAlignment& al) const {
 195     const int32_t absInsertSize = abs(al.InsertSize);
 196     return ( absInsertSize >= MinFragmentLength &&
 197              absInsertSize <= MaxFragmentLength );
 198 }
 199
 200 bool ReadGroupResolver::IsValidOrientation(const BamAlignment& al) const {
 201     const uint16_t currentModelId = CalculateModelType(al) + 1; // convert model type (array index) to ID number
 202     return ( currentModelId == TopModelId || currentModelId == NextTopModelId );
 203 }
 204
 205 void ReadGroupResolver::DetermineTopModels(const string& readGroupName) {
 206
 207     // sort models (from most common to least common)
 208     sort( Models.begin(), Models.end(), std::greater<ModelType>() );
 209
 210     // store top 2 models for later
 211     TopModelId     = Models[0].ID;
 212     NextTopModelId = Models[1].ID;
 213
 214     // make sure that the 2 most common models are some threshold more common
 215     // than the remaining models
 216     const unsigned int activeModelCountSum = Models[0].size() + Models[1].size();
 217     if ( activeModelCountSum == 0 ) return; // skip if no data in this read group
 218     const unsigned int unusedModelCountSum = Models[2].size() + Models[3].size() +
 219                                              Models[4].size() + Models[5].size() +
 220                                              Models[6].size() + Models[7].size();
 221     const double unusedPercentage = (double)unusedModelCountSum / (double)activeModelCountSum;
 222     if ( unusedPercentage > UnusedModelThreshold ) {
 223         cerr << "WARNING: " << readGroupName << " does not have clearly defined 'top models'" << endl
 224              << "         The fraction of alignments in bottom 6 models (" << unusedPercentage
 225              << ") exceeds threshold: " << UnusedModelThreshold << endl;
 226         IsAmbiguous = true;
 227     }
 228
 229     // emit a warning if the best alignment models are non-standard
 230     const bool isModel1Top = (TopModelId == 1) || (NextTopModelId == 1);
 231     const bool isModel2Top = (TopModelId == 2) || (NextTopModelId == 2);
 232     const bool isModel4Top = (TopModelId == 4) || (NextTopModelId == 4);
 233     const bool isModel5Top = (TopModelId == 5) || (NextTopModelId == 5);
 234     const bool isModel6Top = (TopModelId == 6) || (NextTopModelId == 6);
 235     const bool isModel8Top = (TopModelId == 8) || (NextTopModelId == 8);
 236
 237     bool isMatePair  = ( isModel4Top && isModel5Top ? true : false );
 238     bool isPairedEnd = ( isModel2Top && isModel6Top ? true : false );
 239     bool isSolidPair = ( isModel1Top && isModel8Top ? true : false );
 240
 241     if ( !isMatePair && !isPairedEnd && !isSolidPair ) {
 242         cerr << "WARNING: Found a non-standard alignment model configuration. " << endl
 243              << "         Using alignment models " << TopModelId << " & " << NextTopModelId
 244              << endl;
 245     }
 246
 247     // store only the fragments from the best alignment models, then sort
 248     vector<int32_t> fragments;
 249     fragments.reserve( Models[0].size() + Models[1].size() );
 250     fragments.insert( fragments.end(), Models[0].begin(), Models[0].end() );
 251     fragments.insert( fragments.end(), Models[1].begin(), Models[1].end() );
 252     sort ( fragments.begin(), fragments.end() );
 253
 254     // clear out Model fragment data, not needed anymore
 255     Models.clear();
 256
 257     // skip if no fragments found for this read group
 258     if ( fragments.empty() ) {
 259         HasData = false;
 260         return;
 261     } else
 262         HasData = true;
 263
 264     // calculate & store the min,median, & max fragment lengths
 265     const unsigned int numFragmentLengths = fragments.size();
 266     const double halfNonConfidenceInterval = (1.0 - ReadGroupResolver::ConfidenceInterval)/2.0;
 267     const unsigned int minIndex    = (unsigned int)(numFragmentLengths * halfNonConfidenceInterval);
 268     const unsigned int medianIndex = (unsigned int)(numFragmentLengths * 0.5);
 269     const unsigned int maxIndex    = (unsigned int)(numFragmentLengths * (1.0-halfNonConfidenceInterval));
 270
 271     MinFragmentLength    = fragments[minIndex];
 272     MedianFragmentLength = fragments[medianIndex];
 273     MaxFragmentLength    = fragments[maxIndex];
 274 }
 275
 276 void ReadGroupResolver::SetConfidenceInterval(const double& ci) {
 277     ConfidenceInterval = ci;
 278 }
 279
 280 void ReadGroupResolver::SetUnusedModelThreshold(const double& umt) {
 281     UnusedModelThreshold = umt;
 282 }
 283
 284 // --------------------------------------------------------------------------
 285 // ResolveSettings implementation
 286
 287 struct ResolveTool::ResolveSettings {
 288
 289     // flags
 290     bool HasInputBamFile;
 291     bool HasOutputBamFile;
 292     bool IsForceCompression;
 293     bool HasStatsFile;
 294     bool IsMakeStats;
 295     bool IsMarkPairs;
 296     bool IsTwoPass;
 297     bool HasConfidenceInterval;
 298     bool HasUnusedModelThreshold;
 299
 300     // filenames
 301     string InputBamFilename;
 302     string OutputBamFilename;
 303     string StatsFilename;
 304     string ReadNamesFilename; //  ** N.B. - Only used internally, not set from cmdline **
 305
 306     // 'resolve options'
 307     double ConfidenceInterval;
 308     double UnusedModelThreshold;
 309     bool HasForceMarkReadGroups;
 310
 311     // constructor
 312     ResolveSettings(void)
 313         : HasInputBamFile(false)
 314         , HasOutputBamFile(false)
 315         , IsForceCompression(false)
 316         , HasStatsFile(false)
 317         , IsMakeStats(false)
 318         , IsMarkPairs(false)
 319         , IsTwoPass(false)
 320         , HasConfidenceInterval(false)
 321         , HasUnusedModelThreshold(false)
 322         , InputBamFilename(Options::StandardIn())
 323         , OutputBamFilename(Options::StandardOut())
 324         , StatsFilename("")
 325         , ReadNamesFilename(DEFAULT_READNAME_FILE)
 326         , ConfidenceInterval(DEFAULT_CONFIDENCE_INTERVAL)
 327         , UnusedModelThreshold(DEFAULT_UNUSEDMODEL_THRESHOLD)
 328         , HasForceMarkReadGroups(false)
 329     { }
 330 };
 331
 332 // --------------------------------------------------------------------------
 333 // ReadNamesFileReader implementation
 334
 335 struct ResolveTool::ReadNamesFileReader {
 336
 337     // ctor & dtor
 338     ReadNamesFileReader(void) { }
 339     ~ReadNamesFileReader(void) { Close(); }
 340
 341     // main reader interface
 342     public:
 343         void Close(void);
 344         bool Open(const string& filename);
 345         bool Read(map<string, ReadGroupResolver>& readGroups);
 346
 347     // data members
 348     private:
 349         ifstream m_stream;
 350 };
 351
 352 void ResolveTool::ReadNamesFileReader::Close(void) {
 353     if ( m_stream.is_open() )
 354         m_stream.close();
 355 }
 356
 357 bool ResolveTool::ReadNamesFileReader::Open(const string& filename) {
 358
 359     // make sure stream is fresh
 360     Close();
 361
 362     // attempt to open filename, return status
 363     m_stream.open(filename.c_str(), ifstream::in);
 364     return m_stream.good();
 365 }
 366
 367 bool ResolveTool::ReadNamesFileReader::Read(map<string, ReadGroupResolver>& readGroups) {
 368
 369     // up-front sanity check
 370     if ( !m_stream.is_open() ) return false;
 371
 372     // parse read names file
 373     string line;
 374     vector<string> fields;
 375     map<string, ReadGroupResolver>::iterator rgIter;
 376     map<string, ReadGroupResolver>::iterator rgEnd = readGroups.end();
 377     while ( getline(m_stream, line) ) {
 378
 379         // skip if empty line
 380         if ( line.empty() ) continue;
 381
 382         // split line on '\t'
 383         fields = Utilities::Split(line, TAB_CHAR);
 384         if ( fields.size() != 2 ) continue;
 385
 386         // look up resolver for read group
 387         rgIter = readGroups.find( fields[0] );
 388         if ( rgIter == rgEnd ) return false;
 389         ReadGroupResolver& resolver = (*rgIter).second;
 390
 391         // store read name with resolver
 392         resolver.ReadNames.insert( make_pair<string,bool>(fields[1], true) ) ;
 393     }
 394
 395     // if here, return success
 396     return true;
 397 }
 398
 399 // --------------------------------------------------------------------------
 400 // ReadNamesFileWriter implementation
 401
 402 struct ResolveTool::ReadNamesFileWriter {
 403
 404     // ctor & dtor
 405     ReadNamesFileWriter(void) { }
 406     ~ReadNamesFileWriter(void) { Close(); }
 407
 408     // main reader interface
 409     public:
 410         void Close(void);
 411         bool Open(const string& filename);
 412         void Write(const string& readGroupName, const string& readName);
 413
 414     // data members
 415     private:
 416         ofstream m_stream;
 417 };
 418
 419 void ResolveTool::ReadNamesFileWriter::Close(void) {
 420     if ( m_stream.is_open() )
 421         m_stream.close();
 422 }
 423
 424 bool ResolveTool::ReadNamesFileWriter::Open(const string& filename) {
 425
 426     // make sure stream is fresh
 427     Close();
 428
 429     // attempt to open filename, return status
 430     m_stream.open(filename.c_str(), ofstream::out);
 431     return m_stream.good();
 432 }
 433
 434 void ResolveTool::ReadNamesFileWriter::Write(const string& readGroupName,
 435                                              const string& readName)
 436 {
 437     m_stream << readGroupName << TAB_CHAR << readName << endl;
 438 }
 439
 440 // --------------------------------------------------------------------------
 441 // StatsFileReader implementation
 442
 443 struct ResolveTool::StatsFileReader {
 444
 445     // ctor & dtor
 446     public:
 447         StatsFileReader(void) { }
 448         ~StatsFileReader(void) { Close(); }
 449
 450     // main reader interface
 451     public:
 452         void Close(void);
 453         bool Open(const string& filename);
 454         bool Read(ResolveTool::ResolveSettings* settings,
 455                   map<string, ReadGroupResolver>& readGroups);
 456
 457     // internal methods
 458     private:
 459         bool IsComment(const string& line) const;
 460         bool IsWhitespace(const string& line) const;
 461         bool ParseInputLine(const string& line);
 462         bool ParseOptionLine(const string& line, ResolveTool::ResolveSettings* settings);
 463         bool ParseReadGroupLine(const string& line, map<string, ReadGroupResolver>& readGroups);
 464         string SkipCommentsAndWhitespace(void);
 465
 466     // data members
 467     private:
 468         ifstream m_stream;
 469
 470         enum State { None = 0
 471                    , InInput
 472                    , InOptions
 473                    , InReadGroups };
 474 };
 475
 476 void ResolveTool::StatsFileReader::Close(void) {
 477     if ( m_stream.is_open() )
 478         m_stream.close();
 479 }
 480
 481 bool ResolveTool::StatsFileReader::IsComment(const string& line) const {
 482     assert( !line.empty() );
 483     return ( line.at(0) == COMMENT_CHAR );
 484 }
 485
 486 bool ResolveTool::StatsFileReader::IsWhitespace(const string& line) const {
 487     if ( line.empty() )
 488         return true;
 489     return ( isspace(line.at(0)) );
 490 }
 491
 492 bool ResolveTool::StatsFileReader::Open(const string& filename) {
 493
 494     // make sure stream is fresh
 495     Close();
 496
 497     // attempt to open filename, return status
 498     m_stream.open(filename.c_str(), ifstream::in);
 499     return m_stream.good();
 500 }
 501
 502 bool ResolveTool::StatsFileReader::ParseInputLine(const string& /*line*/) {
 503     // input lines are ignored (for now at least), tool will use input from command line
 504     return true;
 505 }
 506
 507 bool ResolveTool::StatsFileReader::ParseOptionLine(const string& line,
 508                                                    ResolveTool::ResolveSettings* settings)
 509 {
 510     // split line into option, value
 511     vector<string> fields = Utilities::Split(line, EQUAL_CHAR);
 512     if ( fields.size() != NUM_OPTIONS_FIELDS )
 513         return false;
 514     const string& option = fields.at(0);
 515     stringstream value(fields.at(1));
 516
 517     // -----------------------------------
 518     // handle option based on keyword
 519
 520     // ConfidenceInterval
 521     if ( option == OPTION_CONFIDENCEINTERVAL ) {
 522         value >> settings->ConfidenceInterval;
 523         settings->HasConfidenceInterval = true;
 524         return true;
 525     }
 526
 527     // UnusedModelThreshold
 528     if ( option == OPTION_UNUSEDMODELTHRESHOLD ) {
 529         value >> settings->UnusedModelThreshold;
 530         settings->HasUnusedModelThreshold = true;
 531         return true;
 532     }
 533
 534     // ForceMarkReadGroups
 535     if ( option == OPTION_FORCEMARKREADGROUPS ) {
 536         value >> settings->HasForceMarkReadGroups;
 537         return true;
 538     }
 539
 540     // otherwise unknown option
 541     cerr << "bamtools resolve ERROR - unrecognized option: " << option << " in stats file" << endl;
 542     return false;
 543 }
 544
 545 bool ResolveTool::StatsFileReader::ParseReadGroupLine(const string& line,
 546                                                       map<string, ReadGroupResolver>& readGroups)
 547 {
 548     // split read group data in to fields
 549     vector<string> fields = Utilities::Split(line, WHITESPACE_CHARS);
 550     if ( fields.size() != NUM_READGROUPS_FIELDS ) return false;
 551
 552     // retrieve RG name
 553     const string& name = fields.at(0);
 554
 555     // populate RG's 'resolver' data
 556     ReadGroupResolver resolver;
 557
 558     stringstream dataStream;
 559     dataStream.str(fields.at(1));
 560     dataStream >> resolver.MedianFragmentLength;
 561     dataStream.clear();
 562
 563     dataStream.str(fields.at(2));
 564     dataStream >> resolver.MinFragmentLength;
 565     dataStream.clear();
 566
 567     dataStream.str(fields.at(3));
 568     dataStream >> resolver.MaxFragmentLength;
 569     dataStream.clear();
 570
 571     dataStream.str(fields.at(4));
 572     dataStream >> resolver.TopModelId;
 573     dataStream.clear();
 574
 575     dataStream.str(fields.at(5));
 576     dataStream >> resolver.NextTopModelId;
 577     dataStream.clear();
 578
 579     resolver.IsAmbiguous = ( fields.at(6) == TRUE_KEYWORD );
 580
 581     // store RG entry and return success
 582     readGroups.insert( make_pair<string, ReadGroupResolver>(name, resolver) );
 583     return true;
 584 }
 585
 586 bool ResolveTool::StatsFileReader::Read(ResolveTool::ResolveSettings* settings,
 587                                         map<string, ReadGroupResolver>& readGroups)
 588 {
 589     // up-front sanity checks
 590     if ( !m_stream.is_open() || settings == 0 )
 591         return false;
 592
 593     // clear out read group data
 594     readGroups.clear();
 595
 596     // initialize state
 597     State currentState = StatsFileReader::None;
 598
 599     // read stats file
 600     string line = SkipCommentsAndWhitespace();
 601     while ( !line.empty() ) {
 602
 603         bool foundError = false;
 604
 605         // switch state on keyword found
 606         if ( Utilities::StartsWith(line, INPUT_TOKEN) )
 607             currentState = StatsFileReader::InInput;
 608         else if ( Utilities::StartsWith(line, OPTIONS_TOKEN) )
 609             currentState = StatsFileReader::InOptions;
 610         else if ( Utilities::StartsWith(line, READGROUPS_TOKEN) )
 611             currentState = StatsFileReader::InReadGroups;
 612
 613         // otherwise parse data line, depending on state
 614         else {
 615             if ( currentState == StatsFileReader::InInput )
 616                 foundError = !ParseInputLine(line);
 617             else if ( currentState == StatsFileReader::InOptions )
 618                 foundError = !ParseOptionLine(line, settings);
 619             else if ( currentState == StatsFileReader::InReadGroups )
 620                 foundError = !ParseReadGroupLine(line, readGroups);
 621             else
 622                 foundError = true;
 623         }
 624
 625         // break out if error found
 626         if ( foundError )
 627             return false;
 628
 629         // get next line
 630         line = SkipCommentsAndWhitespace();
 631     }
 632
 633     // if here, return success
 634     return true;
 635 }
 636
 637 string ResolveTool::StatsFileReader::SkipCommentsAndWhitespace(void) {
 638     string line;
 639     do {
 640         if ( m_stream.eof() )
 641             return string();
 642         getline(m_stream, line);
 643     } while ( IsWhitespace(line) || IsComment(line) );
 644     return line;
 645 }
 646
 647 // --------------------------------------------------------------------------
 648 // StatsFileReader implementation
 649
 650 struct ResolveTool::StatsFileWriter {
 651
 652     // ctor & dtor
 653     public:
 654         StatsFileWriter(void) { }
 655         ~StatsFileWriter(void) { Close(); }
 656
 657     // main reader interface
 658     public:
 659         void Close(void);
 660         bool Open(const string& filename);
 661         bool Write(ResolveTool::ResolveSettings* settings,
 662                    const map<string, ReadGroupResolver>& readGroups);
 663
 664     // internal methods
 665     private:
 666         void WriteHeader(void);
 667         void WriteInput(ResolveTool::ResolveSettings* settings);
 668         void WriteOptions(ResolveTool::ResolveSettings* settings);
 669         void WriteReadGroups(const map<string, ReadGroupResolver>& readGroups);
 670
 671     // data members
 672     private:
 673         ofstream m_stream;
 674 };
 675
 676 void ResolveTool::StatsFileWriter::Close(void) {
 677     if ( m_stream.is_open() )
 678         m_stream.close();
 679 }
 680
 681 bool ResolveTool::StatsFileWriter::Open(const string& filename) {
 682
 683     // make sure stream is fresh
 684     Close();
 685
 686     // attempt to open filename, return status
 687     m_stream.open(filename.c_str(), ofstream::out);
 688     return m_stream.good();
 689 }
 690
 691 bool ResolveTool::StatsFileWriter::Write(ResolveTool::ResolveSettings* settings,
 692                                          const map<string, ReadGroupResolver>& readGroups)
 693 {
 694     // return failure if file not open
 695     if ( !m_stream.is_open() )
 696         return false;
 697
 698     // write stats file elements
 699     WriteHeader();
 700     WriteInput(settings);
 701     WriteOptions(settings);
 702     WriteReadGroups(readGroups);
 703
 704     // return success
 705     return true;
 706 }
 707
 708 void ResolveTool::StatsFileWriter::WriteHeader(void) {
 709
 710     // stringify current bamtools version
 711     stringstream versionStream("");
 712     versionStream << "v"
 713                   << BAMTOOLS_VERSION_MAJOR << "."
 714                   << BAMTOOLS_VERSION_MINOR << "."
 715                   << BAMTOOLS_VERSION_BUILD;
 716
 717     // # bamtools resolve (vX.Y.Z)
 718     // \n
 719
 720     m_stream << COMMENT_CHAR << " bamtools resolve (" << versionStream.str() << ")" << endl
 721              << endl;
 722 }
 723
 724 void ResolveTool::StatsFileWriter::WriteInput(ResolveTool::ResolveSettings* settings) {
 725
 726     // [Input]
 727     // filename
 728     // \n
 729
 730     m_stream << INPUT_TOKEN << endl
 731              << settings->InputBamFilename << endl
 732              << endl;
 733 }
 734
 735 void ResolveTool::StatsFileWriter::WriteOptions(ResolveTool::ResolveSettings* settings) {
 736
 737     // [Options]
 738     // ConfidenceInterval=<double>
 739     // UnusedModelThreshold=<double>
 740     // ForceMarkReadGroups=<true|false>
 741     // \n
 742
 743     m_stream << OPTIONS_TOKEN << endl
 744              << OPTION_CONFIDENCEINTERVAL   << EQUAL_CHAR << settings->ConfidenceInterval << endl
 745              << OPTION_UNUSEDMODELTHRESHOLD << EQUAL_CHAR << settings->UnusedModelThreshold << endl
 746              << OPTION_FORCEMARKREADGROUPS  << EQUAL_CHAR << boolalpha << settings->HasForceMarkReadGroups << endl
 747              << endl;
 748 }
 749
 750 void ResolveTool::StatsFileWriter::WriteReadGroups(const map<string, ReadGroupResolver>& readGroups) {
 751
 752     // [ReadGroups]
 753     // #<name> <medianFL> <minFL> <maxFL> <topModelID> <nextTopModelID> <isAmbiguous?>
 754     m_stream << READGROUPS_TOKEN << endl
 755              << RG_FIELD_DESCRIPTION << endl;
 756
 757     // iterate over read groups
 758     map<string, ReadGroupResolver>::const_iterator rgIter = readGroups.begin();
 759     map<string, ReadGroupResolver>::const_iterator rgEnd  = readGroups.end();
 760     for ( ; rgIter != rgEnd; ++rgIter ) {
 761         const string& name = (*rgIter).first;
 762         const ReadGroupResolver& resolver = (*rgIter).second;
 763
 764         // skip if read group has no data
 765         if ( !resolver.HasData )
 766             continue;
 767
 768         // write read group data
 769         m_stream << name << TAB_CHAR
 770                  << resolver.MedianFragmentLength << TAB_CHAR
 771                  << resolver.MinFragmentLength << TAB_CHAR
 772                  << resolver.MaxFragmentLength << TAB_CHAR
 773                  << resolver.TopModelId << TAB_CHAR
 774                  << resolver.NextTopModelId << TAB_CHAR
 775                  << boolalpha << resolver.IsAmbiguous
 776                  << endl;
 777     }
 778
 779     // extra newline at end
 780     m_stream << endl;
 781 }
 782
 783 // --------------------------------------------------------------------------
 784 // ResolveToolPrivate implementation
 785
 786 struct ResolveTool::ResolveToolPrivate {
 787
 788     // ctor & dtor
 789     public:
 790         ResolveToolPrivate(ResolveTool::ResolveSettings* settings)
 791             : m_settings(settings)
 792         { }
 793         ~ResolveToolPrivate(void) { }
 794
 795     // 'public' interface
 796     public:
 797         bool Run(void);
 798
 799     // internal methods
 800     private:
 801         bool CheckSettings(vector<string>& errors);
 802         bool MakeStats(void);
 803         void ParseHeader(const SamHeader& header);
 804         bool ReadStatsFile(void);
 805         void ResolveAlignment(BamAlignment& al);
 806         bool ResolvePairs(void);
 807         bool WriteStatsFile(void);
 808
 809     // data members
 810     private:
 811         ResolveTool::ResolveSettings* m_settings;
 812         map<string, ReadGroupResolver> m_readGroups;
 813 };
 814
 815 bool ResolveTool::ResolveToolPrivate::CheckSettings(vector<string>& errors) {
 816
 817     // ensure clean slate
 818     errors.clear();
 819
 820     // if MakeStats mode
 821     if ( m_settings->IsMakeStats ) {
 822
 823         // ensure mutex mode
 824         if ( m_settings->IsMarkPairs )
 825             errors.push_back("Cannot run in both -makeStats & -markPairs modes. Please select ONE.");
 826         if ( m_settings->IsTwoPass )
 827             errors.push_back("Cannot run in both -makeStats & -twoPass modes. Please select ONE.");
 828
 829         // error if output BAM options supplied
 830         if ( m_settings->HasOutputBamFile )
 831             errors.push_back("Cannot use -out (output BAM file) in -makeStats mode.");
 832         if ( m_settings->IsForceCompression )
 833             errors.push_back("Cannot use -forceCompression. No output BAM file is being generated.");
 834
 835         // make sure required stats file supplied
 836         if ( !m_settings->HasStatsFile )
 837             errors.push_back("Ouptut stats filename required for -makeStats mode. Please specify one using -stats option.");
 838
 839         // check for UseStats options
 840         if ( m_settings->HasForceMarkReadGroups )
 841             errors.push_back("Cannot use -forceMarkReadGroups. -markPairs options are DISABLED in -makeStats mode.");
 842     }
 843
 844     // if MarkPairs mode
 845     else if ( m_settings->IsMarkPairs ) {
 846
 847         // ensure mutex mode
 848         if ( m_settings->IsMakeStats )
 849             errors.push_back("Cannot run in both -makeStats & -markPairs modes. Please select ONE.");
 850         if ( m_settings->IsTwoPass )
 851             errors.push_back("Cannot run in both -markPairs & -twoPass modes. Please select ONE.");
 852
 853         // make sure required stats file supplied
 854         if ( !m_settings->HasStatsFile )
 855             errors.push_back("Input stats filename required for -markPairs mode. Please specify one using -stats option.");
 856
 857         // check for MakeStats options
 858         if ( m_settings->HasConfidenceInterval )
 859             errors.push_back("Cannot use -ci. -makeStats options are DISABLED is -markPairs mode.");
 860     }
 861
 862     // if TwoPass mode
 863     else if ( m_settings->IsTwoPass ) {
 864
 865         // ensure mutex mode
 866         if ( m_settings->IsMakeStats )
 867             errors.push_back("Cannot run in both -makeStats & -twoPass modes. Please select ONE.");
 868         if ( m_settings->IsMarkPairs )
 869             errors.push_back("Cannot run in both -markPairs & -twoPass modes. Please select ONE.");
 870
 871         // make sure input is file not stdin
 872         if ( !m_settings->HasInputBamFile || m_settings->InputBamFilename == Options::StandardIn() )
 873             errors.push_back("Cannot run -twoPass mode with BAM data from stdin. Please specify existing file using -in option.");
 874     }
 875
 876     // no mode selected
 877     else
 878         errors.push_back("No resolve mode specified. Please select ONE of the following: -makeStats, -markPairs, or -twoPass. See help for more info.");
 879
 880     // boundary checks on values
 881     if ( m_settings->HasConfidenceInterval ) {
 882         if ( m_settings->ConfidenceInterval < 0.0 || m_settings->ConfidenceInterval > 1.0 )
 883             errors.push_back("Invalid confidence interval. Must be between 0 and 1");
 884     }
 885     if ( m_settings->HasUnusedModelThreshold ) {
 886         if ( m_settings->UnusedModelThreshold < 0.0 || m_settings->UnusedModelThreshold > 1.0 )
 887             errors.push_back("Invalid unused model threshold. Must be between 0 and 1");
 888     }
 889
 890     // return success if no errors found
 891     return ( errors.empty() );
 892 }
 893
 894 bool ResolveTool::ResolveToolPrivate::MakeStats(void) {
 895
 896     // pull resolver settings from command-line settings
 897     ReadGroupResolver::SetConfidenceInterval(m_settings->ConfidenceInterval);
 898     ReadGroupResolver::SetUnusedModelThreshold(m_settings->UnusedModelThreshold);
 899
 900     // open our BAM reader
 901     BamReader bamReader;
 902     if ( !bamReader.Open(m_settings->InputBamFilename) ) {
 903         cerr << "bamtools resolve ERROR: could not open input BAM file: "
 904              << m_settings->InputBamFilename << endl;
 905         return false;
 906     }
 907
 908     // retrieve header & parse for read groups
 909     const SamHeader& header = bamReader.GetHeader();
 910     ParseHeader(header);
 911
 912     // open ReadNamesFileWriter
 913     ResolveTool::ReadNamesFileWriter readNamesWriter;
 914     if ( !readNamesWriter.Open(m_settings->ReadNamesFilename) ) {
 915         cerr << "bamtools resolve ERROR: could not open (temp) output read names file: "
 916              << m_settings->ReadNamesFilename << endl;
 917         bamReader.Close();
 918         return false;
 919     }
 920
 921     // read through BAM file
 922     BamAlignment al;
 923     string readGroup("");
 924     map<string, ReadGroupResolver>::iterator rgIter;
 925     map<string, bool>::iterator readNameIter;
 926     while ( bamReader.GetNextAlignmentCore(al) ) {
 927
 928         // skip if alignment is not paired, mapped, nor mate is mapped
 929         if ( !al.IsPaired() || !al.IsMapped() || !al.IsMateMapped() )
 930             continue;
 931
 932         // skip if alignment & mate not on same reference sequence
 933         if ( al.RefID != al.MateRefID ) continue;
 934
 935         // flesh out the char data, so we can retrieve its read group ID
 936         al.BuildCharData();
 937
 938         // get read group from alignment (OK if empty)
 939         readGroup.clear();
 940         al.GetTag(READ_GROUP_TAG, readGroup);
 941
 942         // look up resolver for read group
 943         rgIter = m_readGroups.find(readGroup);
 944         if ( rgIter == m_readGroups.end() )  {
 945             cerr << "bamtools resolve ERROR - unable to calculate stats, unknown read group encountered: "
 946                  << readGroup << endl;
 947             bamReader.Close();
 948             return false;
 949         }
 950         ReadGroupResolver& resolver = (*rgIter).second;
 951
 952         // determine unique-ness of current alignment
 953         const bool isCurrentMateUnique = ( al.MapQuality != 0 );
 954
 955         // look up read name
 956         readNameIter = resolver.ReadNames.find(al.Name);
 957
 958         // if read name found (current alignment's mate already parsed)
 959         if ( readNameIter != resolver.ReadNames.end() ) {
 960
 961             // if both unique mates are unique, store read name & insert size for later
 962             const bool isStoredMateUnique  = (*readNameIter).second;
 963             if ( isCurrentMateUnique && isStoredMateUnique ) {
 964
 965                 // save read name in temp file as candidates for later pair marking
 966                 readNamesWriter.Write(readGroup, al.Name);
 967
 968                 // determine model type & store fragment length for stats calculation
 969                 const uint16_t currentModelType = CalculateModelType(al);
 970                 assert( currentModelType != ModelType::DUMMY_ID );
 971                 resolver.Models[currentModelType].push_back( abs(al.InsertSize) );
 972             }
 973
 974             // unique or not, remove read name from map
 975             resolver.ReadNames.erase(readNameIter);
 976         }
 977
 978         // if read name not found, store new entry
 979         else resolver.ReadNames.insert( make_pair<string, bool>(al.Name, isCurrentMateUnique) );
 980     }
 981
 982     // close files
 983     readNamesWriter.Close();
 984     bamReader.Close();
 985
 986     // iterate back through read groups
 987     map<string, ReadGroupResolver>::iterator rgEnd  = m_readGroups.end();
 988     for ( rgIter = m_readGroups.begin(); rgIter != rgEnd; ++rgIter ) {
 989         const string& name = (*rgIter).first;
 990         ReadGroupResolver& resolver = (*rgIter).second;
 991
 992         // calculate acceptable orientation & insert sizes for this read group
 993         resolver.DetermineTopModels(name);
 994
 995         // clear out left over read names
 996         // (these have mates that did not pass filters or were already removed as non-unique)
 997         resolver.ReadNames.clear();
 998     }
 999
1000     // if we get here, return success
1001     return true;
1002 }
1003
1004 void ResolveTool::ResolveToolPrivate::ParseHeader(const SamHeader& header) {
1005
1006     // iterate over header read groups, creating a 'resolver' for each
1007     SamReadGroupConstIterator rgIter = header.ReadGroups.ConstBegin();
1008     SamReadGroupConstIterator rgEnd  = header.ReadGroups.ConstEnd();
1009     for ( ; rgIter != rgEnd; ++rgIter ) {
1010         const SamReadGroup& rg = (*rgIter);
1011         m_readGroups.insert( make_pair<string, ReadGroupResolver>(rg.ID, ReadGroupResolver()) );
1012     }
1013 }
1014
1015 bool ResolveTool::ResolveToolPrivate::ReadStatsFile(void) {
1016
1017     // skip if no filename provided
1018     if ( m_settings->StatsFilename.empty() )
1019         return false;
1020
1021     // attempt to open stats file
1022     ResolveTool::StatsFileReader statsReader;
1023     if ( !statsReader.Open(m_settings->StatsFilename) ) {
1024         cerr << "bamtools resolve ERROR - could not open stats file: "
1025              << m_settings->StatsFilename << " for reading" << endl;
1026         return false;
1027     }
1028
1029     // attempt to read stats data
1030     if ( !statsReader.Read(m_settings, m_readGroups) ) {
1031         cerr << "bamtools resolve ERROR - could not parse stats file: "
1032              << m_settings->StatsFilename << " for data" << endl;
1033         return false;
1034     }
1035
1036     // return success
1037     return true;
1038 }
1039
1040 void ResolveTool::ResolveToolPrivate::ResolveAlignment(BamAlignment& al) {
1041
1042     // clear proper-pair flag
1043     al.SetIsProperPair(false);
1044
1045     // quit check if alignment is not from paired-end read
1046     if ( !al.IsPaired() ) return;
1047
1048     // quit check if either alignment or its mate are unmapped
1049     if ( !al.IsMapped() || !al.IsMateMapped() ) return;
1050
1051     // quit check if alignment & its mate are on differenct references
1052     if ( al.RefID != al.MateRefID ) return;
1053
1054     // quit check if map quality is 0
1055     if ( al.MapQuality == 0 ) return;
1056
1057     // get read group from alignment
1058     // empty string if not found, this is OK - we handle empty read group case
1059     string readGroupName("");
1060     al.GetTag(READ_GROUP_TAG, readGroupName);
1061
1062     // look up read group's 'resolver'
1063     map<string, ReadGroupResolver>::iterator rgIter = m_readGroups.find(readGroupName);
1064     if ( rgIter == m_readGroups.end() ) {
1065         cerr << "bamtools resolve ERROR - read group found that was not in header: "
1066              << readGroupName << endl;
1067         exit(1);
1068     }
1069     const ReadGroupResolver& resolver = (*rgIter).second;
1070
1071     // quit check if pairs are not in proper orientation (can differ for each RG)
1072     if ( !resolver.IsValidOrientation(al) ) return;
1073
1074     // quit check if pairs are not within "reasonable" distance (can differ for each RG)
1075     if ( !resolver.IsValidInsertSize(al) ) return;
1076
1077     // quit check if alignment is not a "candidate proper pair"
1078     map<string, bool>::const_iterator readNameIter;
1079     readNameIter = resolver.ReadNames.find(al.Name);
1080     if ( readNameIter == resolver.ReadNames.end() )
1081         return;
1082
1083     // if we get here, alignment is OK - set 'proper pair' flag
1084     al.SetIsProperPair(true);
1085 }
1086
1087 bool ResolveTool::ResolveToolPrivate::ResolvePairs(void) {
1088
1089     // open file containing read names of candidate proper pairs
1090     ResolveTool::ReadNamesFileReader readNamesReader;
1091     if ( !readNamesReader.Open(m_settings->ReadNamesFilename) ) {
1092         cerr << "bamtools resolve ERROR: could not open (temp) inputput read names file: "
1093              << m_settings->ReadNamesFilename << endl;
1094         return false;
1095     }
1096
1097     // parse read names (matching with corresponding read groups)
1098     if ( !readNamesReader.Read(m_readGroups) ) {
1099         cerr << "bamtools resolve ERROR: could not read candidate read names from file: "
1100              << m_settings->ReadNamesFilename << endl;
1101         readNamesReader.Close();
1102         return false;
1103     }
1104
1105     // close read name file reader & delete temp file
1106     readNamesReader.Close();
1107     if ( remove(m_settings->ReadNamesFilename.c_str()) != 0 ) {
1108         cerr << "bamtools resolve WARNING: could not delete temp file: "
1109              << m_settings->ReadNamesFilename << endl;
1110     }
1111
1112     // open our BAM reader
1113     BamReader reader;
1114     if ( !reader.Open(m_settings->InputBamFilename) ) {
1115         cerr << "bamtools resolve ERROR: could not open input BAM file: "
1116              << m_settings->InputBamFilename << endl;
1117         return false;
1118     }
1119
1120     // retrieve header & reference dictionary info
1121     const SamHeader& header = reader.GetHeader();
1122     const RefVector& references = reader.GetReferenceData();
1123
1124     // determine compression mode for BamWriter
1125     bool writeUncompressed = ( m_settings->OutputBamFilename == Options::StandardOut() &&
1126                                !m_settings->IsForceCompression );
1127     BamWriter::CompressionMode compressionMode = BamWriter::Compressed;
1128     if ( writeUncompressed ) compressionMode = BamWriter::Uncompressed;
1129
1130     // open BamWriter
1131     BamWriter writer;
1132     writer.SetCompressionMode(compressionMode);
1133     if ( !writer.Open(m_settings->OutputBamFilename, header, references) ) {
1134         cerr << "bamtools resolve ERROR: could not open "
1135              << m_settings->OutputBamFilename << " for writing." << endl;
1136         reader.Close();
1137         return false;
1138     }
1139
1140     // plow through alignments, setting/clearing 'proper pair' flag
1141     // and writing to new output BAM file
1142     BamAlignment al;
1143     while ( reader.GetNextAlignment(al) ) {
1144         ResolveAlignment(al);
1145         writer.SaveAlignment(al);
1146     }
1147
1148     // clean up & return success
1149     reader.Close();
1150     writer.Close();
1151     return true;
1152 }
1153
1154 bool ResolveTool::ResolveToolPrivate::Run(void) {
1155
1156     // verify that command line settings are acceptable
1157     vector<string> errors;
1158     if ( !CheckSettings(errors) ) {
1159         cerr << "bamtools resolve ERROR - invalid settings: " << endl;
1160         vector<string>::const_iterator errorIter = errors.begin();
1161         vector<string>::const_iterator errorEnd  = errors.end();
1162         for ( ; errorIter != errorEnd; ++errorIter )
1163             cerr << (*errorIter) << endl;
1164         return false;
1165     }
1166
1167     // initialize read group map with default (empty name) read group
1168     m_readGroups.insert( make_pair<string, ReadGroupResolver>("", ReadGroupResolver()) );
1169
1170     // init readname filename
1171     // uses (adjusted) stats filename if provided (req'd for makeStats, markPairs modes; optional for twoPass)
1172     // else keep default filename
1173     if ( m_settings->HasStatsFile )
1174         m_settings->ReadNamesFilename = m_settings->StatsFilename + READNAME_FILE_SUFFIX;
1175
1176     // -makeStats mode
1177     if ( m_settings->IsMakeStats ) {
1178
1179         // generate stats data
1180         if ( !MakeStats() ) {
1181             cerr << "bamtools resolve ERROR - could not generate stats" << endl;
1182             return false;
1183         }
1184
1185         // write stats to file
1186         if ( !WriteStatsFile() ) {
1187             cerr << "bamtools resolve ERROR - could not write stats file: "
1188                  << m_settings->StatsFilename << endl;
1189             return false;
1190         }
1191     }
1192
1193     // -markPairs mode
1194     else if ( m_settings->IsMarkPairs ) {
1195
1196         // read stats from file
1197         if ( !ReadStatsFile() ) {
1198             cerr << "bamtools resolve ERROR - could not read stats file: "
1199                  << m_settings->StatsFilename << endl;
1200             return false;
1201         }
1202
1203         // do paired-end resolution
1204         if ( !ResolvePairs() ) {
1205             cerr << "bamtools resolve ERROR - could not resolve pairs" << endl;
1206             return false;
1207         }
1208     }
1209
1210     // -twoPass mode
1211     else {
1212
1213         // generate stats data
1214         if ( !MakeStats() ) {
1215             cerr << "bamtools resolve ERROR - could not generate stats" << endl;
1216             return false;
1217         }
1218
1219         // if stats file requested
1220         if ( m_settings->HasStatsFile ) {
1221
1222             // write stats to file
1223             // emit warning if write fails, but paired-end resolution should be allowed to proceed
1224             if ( !WriteStatsFile() )
1225                 cerr << "bamtools resolve WARNING - could not write stats file: "
1226                      << m_settings->StatsFilename << endl;
1227         }
1228
1229         // do paired-end resolution
1230         if ( !ResolvePairs() ) {
1231             cerr << "bamtools resolve ERROR - could not resolve pairs" << endl;
1232             return false;
1233         }
1234     }
1235
1236     // return success
1237     return true;
1238 }
1239
1240 bool ResolveTool::ResolveToolPrivate::WriteStatsFile(void) {
1241
1242     // skip if no filename provided
1243     if ( m_settings->StatsFilename.empty() )
1244         return false;
1245
1246     // attempt to open stats file
1247     ResolveTool::StatsFileWriter statsWriter;
1248     if ( !statsWriter.Open(m_settings->StatsFilename) ) {
1249         cerr << "bamtools resolve ERROR - could not open stats file: "
1250              << m_settings->StatsFilename << " for writing" << endl;
1251         return false;
1252     }
1253
1254     // attempt to write stats data
1255     if ( !statsWriter.Write(m_settings, m_readGroups) ) {
1256         cerr << "bamtools resolve ERROR - could not write stats file: "
1257              << m_settings->StatsFilename << " for data" << endl;
1258         return false;
1259     }
1260
1261     // return success
1262     return true;
1263 }
1264
1265 // --------------------------------------------------------------------------
1266 // ResolveTool implementation
1267
1268 ResolveTool::ResolveTool(void)
1269     : AbstractTool()
1270     , m_settings(new ResolveSettings)
1271     , m_impl(0)
1272 {
1273     // set description texts
1274     const string programDescription = "resolves paired-end reads (marking the IsProperPair flag as needed)";
1275     const string programUsage = "<mode> [options] [-in <filename>] [-out <filename> | [-forceCompression] ] [-stats <filename>]";
1276     const string inputBamDescription = "the input BAM file(s)";
1277     const string outputBamDescription = "the output BAM file";
1278     const string statsFileDescription = "input/output stats file, depending on selected mode (see below). "
1279             "This file is human-readable, storing fragment length data generated per read group, as well as "
1280             "the options used to configure the -makeStats mode";
1281     const string forceCompressionDescription = "if results are sent to stdout (like when piping to another tool), "
1282             "default behavior is to leave output uncompressed."
1283             "Use this flag to override and force compression. This feature is disabled in -makeStats mode.";
1284     const string makeStatsDescription = "generates a fragment-length stats file from the input BAM. "
1285             "Data is written to file specified using the -stats option. "
1286             "MarkPairs Mode Settings are DISABLED.";
1287     const string markPairsDescription = "generates an output BAM with alignments marked with proper-pair status. "
1288             "Stats data is read from file specified using the -stats option. "
1289             "MakeStats Mode Settings are DISABLED";
1290     const string twoPassDescription = "combines the -makeStats & -markPairs modes into a single command. "
1291             "However, due to the two-pass nature of paired-end resolution, piping BAM data via stdin is DISABLED. "
1292             "You must supply an explicit input BAM file. Output BAM may be piped to stdout, however, if desired. "
1293             "All MakeStats & MarkPairs Mode Settings are available. "
1294             "The intermediate stats file is not necessary, but if the -stats options is used, then one will be generated. "
1295             "You may find this useful for documentation purposes.";
1296     const string confidenceIntervalDescription = "confidence interval. Set min/max fragment lengths such that we capture "
1297             "this fraction of pairs";
1298     const string unusedModelThresholdDescription = "unused model threshold. The resolve tool considers 8 possible orientation models "
1299             "for pairs. The top 2 are selected for later use when actually marking alignments. This value determines the "
1300             "cutoff for marking a read group as ambiguous. Meaning that if the ratio of the number of alignments from bottom 6 models "
1301             "to the top 2 is greater than this threshold, then the read group is flagged as ambiguous. By default, NO alignments "
1302             "from ambiguous read groups will be marked as proper pairs. You may override this behavior with the -force option "
1303             "in -markPairs mode";
1304     const string forceMarkDescription = "forces all read groups to be marked according to their top 2 'orientation models'. "
1305             "When generating stats, the 2 (out of 8 possible) models with the most observations are chosen as the top models for each read group. "
1306             "If the remaining 6 models account for more than some threshold ([default=10%], see -umt), then the read group is marked as ambiguous. "
1307             "The default behavior is that for an ambiguous read group, NONE of its alignments are marked as proper-pairs. "
1308             "By setting this option, a read group's ambiguity flag will be ignored, and all of its alignments will be compared to the top 2 models.";
1309
1310     // set program details
1311     Options::SetProgramInfo("bamtools resolve", programDescription, programUsage);
1312
1313     // set up I/O options
1314     OptionGroup* IO_Opts = Options::CreateOptionGroup("Input & Output");
1315     Options::AddValueOption("-in",  "BAM filename", inputBamDescription, "",
1316                             m_settings->HasInputBamFile, m_settings->InputBamFilename,
1317                             IO_Opts, Options::StandardIn());
1318     Options::AddValueOption("-out", "BAM filename", outputBamDescription, "",
1319                             m_settings->HasOutputBamFile, m_settings->OutputBamFilename,
1320                             IO_Opts, Options::StandardOut());
1321     Options::AddValueOption("-stats", "STATS filename", statsFileDescription, "",
1322                             m_settings->HasStatsFile, m_settings->StatsFilename, IO_Opts);
1323     Options::AddOption("-forceCompression", forceCompressionDescription,
1324                        m_settings->IsForceCompression, IO_Opts);
1325
1326     OptionGroup* ModeOpts = Options::CreateOptionGroup("Resolve Modes (must select ONE of the following)");
1327     Options::AddOption("-makeStats", makeStatsDescription, m_settings->IsMakeStats, ModeOpts);
1328     Options::AddOption("-markPairs", markPairsDescription, m_settings->IsMarkPairs, ModeOpts);
1329     Options::AddOption("-twoPass",   twoPassDescription,   m_settings->IsTwoPass,   ModeOpts);
1330
1331     OptionGroup* MakeStatsOpts = Options::CreateOptionGroup("MakeStats Mode Options (disabled in -markPairs mode)");
1332     Options::AddValueOption("-ci", "double", confidenceIntervalDescription, "",
1333                             m_settings->HasConfidenceInterval, m_settings->ConfidenceInterval, MakeStatsOpts);
1334     Options::AddValueOption("-umt", "double", unusedModelThresholdDescription, "",
1335                             m_settings->HasUnusedModelThreshold, m_settings->UnusedModelThreshold, MakeStatsOpts);
1336
1337     OptionGroup* MarkPairsOpts = Options::CreateOptionGroup("MarkPairs Mode Options (disabled in -makeStats mode)");
1338     Options::AddOption("-force", forceMarkDescription, m_settings->HasForceMarkReadGroups, MarkPairsOpts);
1339 }
1340
1341 ResolveTool::~ResolveTool(void) {
1342
1343     delete m_settings;
1344     m_settings = 0;
1345
1346     delete m_impl;
1347     m_impl = 0;
1348 }
1349
1350 int ResolveTool::Help(void) {
1351     Options::DisplayHelp();
1352     return 0;
1353 }
1354
1355 int ResolveTool::Run(int argc, char* argv[]) {
1356
1357     // parse command line arguments
1358     Options::Parse(argc, argv, 1);
1359
1360     // initialize ResolveTool
1361     m_impl = new ResolveToolPrivate(m_settings);
1362
1363     // run ResolveTool, return success/failure
1364     if ( m_impl->Run() )
1365         return 0;
1366     else
1367         return 1;
1368 }