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