src/toolkit/bamtools_convert.cpp

   1 // ***************************************************************************
   2 // bamtools_convert.cpp (c) 2010 Derek Barnett, Erik Garrison
   3 // Marth Lab, Department of Biology, Boston College
   4 // ---------------------------------------------------------------------------
   5 // Last modified: 10 December 2012
   6 // ---------------------------------------------------------------------------
   7 // Converts between BAM and a number of other formats
   8 // ***************************************************************************
   9
  10 #include "bamtools_convert.h"
  11
  12 #include <api/BamConstants.h>
  13 #include <api/BamMultiReader.h>
  14 #include <utils/bamtools_fasta.h>
  15 #include <utils/bamtools_options.h>
  16 #include <utils/bamtools_pileup_engine.h>
  17 #include <utils/bamtools_utilities.h>
  18 using namespace BamTools;
  19
  20 #include <fstream>
  21 #include <iostream>
  22 #include <sstream>
  23 #include <string>
  24 #include <vector>
  25 using namespace std;
  26
  27 namespace BamTools {
  28
  29 // ---------------------------------------------
  30 // ConvertTool constants
  31
  32 // supported conversion format command-line names
  33 static const string FORMAT_BED    = "bed";
  34 static const string FORMAT_FASTA  = "fasta";
  35 static const string FORMAT_FASTQ  = "fastq";
  36 static const string FORMAT_JSON   = "json";
  37 static const string FORMAT_SAM    = "sam";
  38 static const string FORMAT_PILEUP = "pileup";
  39 static const string FORMAT_YAML   = "yaml";
  40
  41 // other constants
  42 static const unsigned int FASTA_LINE_MAX = 50;
  43
  44 // ---------------------------------------------
  45 // ConvertPileupFormatVisitor declaration
  46
  47 class ConvertPileupFormatVisitor : public PileupVisitor {
  48
  49     // ctor & dtor
  50     public:
  51         ConvertPileupFormatVisitor(const RefVector& references,
  52                                    const string& fastaFilename,
  53                                    const bool isPrintingMapQualities,
  54                                    ostream* out);
  55         ~ConvertPileupFormatVisitor(void);
  56
  57     // PileupVisitor interface implementation
  58     public:
  59         void Visit(const PileupPosition& pileupData);
  60
  61     // data members
  62     private:
  63         Fasta     m_fasta;
  64         bool      m_hasFasta;
  65         bool      m_isPrintingMapQualities;
  66         ostream*  m_out;
  67         RefVector m_references;
  68 };
  69
  70 } // namespace BamTools
  71
  72 // ---------------------------------------------
  73 // ConvertSettings implementation
  74
  75 struct ConvertTool::ConvertSettings {
  76
  77     // flag
  78     bool HasInput;
  79     bool HasInputFilelist;
  80     bool HasOutput;
  81     bool HasFormat;
  82     bool HasRegion;
  83
  84     // pileup flags
  85     bool HasFastaFilename;
  86     bool IsOmittingSamHeader;
  87     bool IsPrintingPileupMapQualities;
  88
  89     // options
  90     vector<string> InputFiles;
  91     string InputFilelist;
  92     string OutputFilename;
  93     string Format;
  94     string Region;
  95
  96     // pileup options
  97     string FastaFilename;
  98
  99     // constructor
 100     ConvertSettings(void)
 101         : HasInput(false)
 102         , HasInputFilelist(false)
 103         , HasOutput(false)
 104         , HasFormat(false)
 105         , HasRegion(false)
 106         , HasFastaFilename(false)
 107         , IsOmittingSamHeader(false)
 108         , IsPrintingPileupMapQualities(false)
 109         , OutputFilename(Options::StandardOut())
 110         , FastaFilename("")
 111     { }
 112 };
 113
 114 // ---------------------------------------------
 115 // ConvertToolPrivate implementation
 116
 117 struct ConvertTool::ConvertToolPrivate {
 118
 119     // ctor & dtor
 120     public:
 121         ConvertToolPrivate(ConvertTool::ConvertSettings* settings)
 122             : m_settings(settings)
 123             , m_out(cout.rdbuf())
 124         { }
 125
 126         ~ConvertToolPrivate(void) { }
 127
 128     // interface
 129     public:
 130         bool Run(void);
 131
 132     // internal methods
 133     private:
 134         void PrintBed(const BamAlignment& a);
 135         void PrintFasta(const BamAlignment& a);
 136         void PrintFastq(const BamAlignment& a);
 137         void PrintJson(const BamAlignment& a);
 138         void PrintSam(const BamAlignment& a);
 139         void PrintYaml(const BamAlignment& a);
 140
 141         // special case - uses the PileupEngine
 142         bool RunPileupConversion(BamMultiReader* reader);
 143
 144     // data members
 145     private:
 146         ConvertTool::ConvertSettings* m_settings;
 147         RefVector m_references;
 148         ostream m_out;
 149 };
 150
 151 bool ConvertTool::ConvertToolPrivate::Run(void) {
 152
 153     // ------------------------------------
 154     // initialize conversion input/output
 155
 156     // set to default input if none provided
 157     if ( !m_settings->HasInput && !m_settings->HasInputFilelist )
 158         m_settings->InputFiles.push_back(Options::StandardIn());
 159
 160     // add files in the filelist to the input file list
 161     if ( m_settings->HasInputFilelist ) {
 162
 163         ifstream filelist(m_settings->InputFilelist.c_str(), ios::in);
 164         if ( !filelist.is_open() ) {
 165             cerr << "bamtools convert ERROR: could not open input BAM file list... Aborting." << endl;
 166             return false;
 167         }
 168
 169         string line;
 170         while ( getline(filelist, line) )
 171             m_settings->InputFiles.push_back(line);
 172     }
 173
 174     // open input files
 175     BamMultiReader reader;
 176     if ( !reader.Open(m_settings->InputFiles) ) {
 177         cerr << "bamtools convert ERROR: could not open input BAM file(s)... Aborting." << endl;
 178         return false;
 179     }
 180
 181     // if input is not stdin & a region is provided, look for index files
 182     if ( m_settings->HasInput && m_settings->HasRegion ) {
 183         if ( !reader.LocateIndexes() ) {
 184             cerr << "bamtools convert ERROR: could not locate index file(s)... Aborting." << endl;
 185             return false;
 186         }
 187     }
 188
 189     // retrieve reference data
 190     m_references = reader.GetReferenceData();
 191
 192     // set region if specified
 193     BamRegion region;
 194     if ( m_settings->HasRegion ) {
 195         if ( Utilities::ParseRegionString(m_settings->Region, reader, region) ) {
 196
 197             if ( reader.HasIndexes() ) {
 198                 if ( !reader.SetRegion(region) ) {
 199                     cerr << "bamtools convert ERROR: set region failed. Check that REGION describes a valid range" << endl;
 200                     reader.Close();
 201                     return false;
 202                 }
 203             }
 204
 205         } else {
 206             cerr << "bamtools convert ERROR: could not parse REGION: " << m_settings->Region << endl;
 207             cerr << "Check that REGION is in valid format (see documentation) and that the coordinates are valid"
 208                  << endl;
 209             reader.Close();
 210             return false;
 211         }
 212     }
 213
 214     // if output file given
 215     ofstream outFile;
 216     if ( m_settings->HasOutput ) {
 217
 218         // open output file stream
 219         outFile.open(m_settings->OutputFilename.c_str());
 220         if ( !outFile ) {
 221             cerr << "bamtools convert ERROR: could not open " << m_settings->OutputFilename
 222                  << " for output" << endl;
 223             return false;
 224         }
 225
 226         // set m_out to file's streambuf
 227         m_out.rdbuf(outFile.rdbuf());
 228     }
 229
 230     // -------------------------------------
 231     // do conversion based on format
 232
 233      bool convertedOk = true;
 234
 235     // pileup is special case
 236     // conversion not done per alignment, like the other formats
 237     if ( m_settings->Format == FORMAT_PILEUP )
 238         convertedOk = RunPileupConversion(&reader);
 239
 240     // all other formats
 241     else {
 242
 243         bool formatError = false;
 244
 245         // set function pointer to proper conversion method
 246         void (BamTools::ConvertTool::ConvertToolPrivate::*pFunction)(const BamAlignment&) = 0;
 247         if      ( m_settings->Format == FORMAT_BED )   pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintBed;
 248         else if ( m_settings->Format == FORMAT_FASTA ) pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintFasta;
 249         else if ( m_settings->Format == FORMAT_FASTQ ) pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintFastq;
 250         else if ( m_settings->Format == FORMAT_JSON )  pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintJson;
 251         else if ( m_settings->Format == FORMAT_SAM )   pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintSam;
 252         else if ( m_settings->Format == FORMAT_YAML )  pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintYaml;
 253         else {
 254             cerr << "bamtools convert ERROR: unrecognized format: " << m_settings->Format << endl;
 255             cerr << "Please see documentation for list of supported formats " << endl;
 256             formatError = true;
 257             convertedOk = false;
 258         }
 259
 260         // if format selected ok
 261         if ( !formatError ) {
 262
 263             // if SAM format & not omitting header, print SAM header first
 264             if ( (m_settings->Format == FORMAT_SAM) && !m_settings->IsOmittingSamHeader )
 265                 m_out << reader.GetHeaderText();
 266
 267             // iterate through file, doing conversion
 268             BamAlignment a;
 269             while ( reader.GetNextAlignment(a) )
 270                 (this->*pFunction)(a);
 271
 272             // set flag for successful conversion
 273             convertedOk = true;
 274         }
 275     }
 276
 277     // ------------------------
 278     // clean up & exit
 279     reader.Close();
 280     if ( m_settings->HasOutput )
 281         outFile.close();
 282     return convertedOk;
 283 }
 284
 285 // ----------------------------------------------------------
 286 // Conversion/output methods
 287 // ----------------------------------------------------------
 288
 289 void ConvertTool::ConvertToolPrivate::PrintBed(const BamAlignment& a) {
 290
 291     // tab-delimited, 0-based half-open
 292     // (e.g. a 50-base read aligned to pos 10 could have BED coordinates (10, 60) instead of BAM coordinates (10, 59) )
 293     // <chromName> <chromStart> <chromEnd> <readName> <score> <strand>
 294
 295     m_out << m_references.at(a.RefID).RefName << "\t"
 296           << a.Position << "\t"
 297           << a.GetEndPosition() << "\t"
 298           << a.Name << "\t"
 299           << a.MapQuality << "\t"
 300           << (a.IsReverseStrand() ? "-" : "+") << endl;
 301 }
 302
 303 // print BamAlignment in FASTA format
 304 // N.B. - uses QueryBases NOT AlignedBases
 305 void ConvertTool::ConvertToolPrivate::PrintFasta(const BamAlignment& a) {
 306
 307     // >BamAlignment.Name
 308     // BamAlignment.QueryBases (up to FASTA_LINE_MAX bases per line)
 309     // ...
 310     //
 311     // N.B. - QueryBases are reverse-complemented if aligned to reverse strand
 312
 313     // print header
 314     m_out << ">" << a.Name << endl;
 315
 316     // handle reverse strand alignment - bases
 317     string sequence = a.QueryBases;
 318     if ( a.IsReverseStrand() )
 319         Utilities::ReverseComplement(sequence);
 320
 321     // if sequence fits on single line
 322     if ( sequence.length() <= FASTA_LINE_MAX )
 323         m_out << sequence << endl;
 324
 325     // else split over multiple lines
 326     else {
 327
 328         size_t position = 0;
 329         size_t seqLength = sequence.length(); // handle reverse strand alignment - bases & qualitiesth();
 330
 331         // write subsequences to each line
 332         while ( position < (seqLength - FASTA_LINE_MAX) ) {
 333             m_out << sequence.substr(position, FASTA_LINE_MAX) << endl;
 334             position += FASTA_LINE_MAX;
 335         }
 336
 337         // write final subsequence
 338         m_out << sequence.substr(position) << endl;
 339     }
 340 }
 341
 342 // print BamAlignment in FASTQ format
 343 // N.B. - uses QueryBases NOT AlignedBases
 344 void ConvertTool::ConvertToolPrivate::PrintFastq(const BamAlignment& a) {
 345
 346     // @BamAlignment.Name
 347     // BamAlignment.QueryBases
 348     // +
 349     // BamAlignment.Qualities
 350     //
 351     // N.B. - QueryBases are reverse-complemented (& Qualities reversed) if aligned to reverse strand .
 352     //        Name is appended "/1" or "/2" if paired-end, to reflect which mate this entry is.
 353
 354     // handle paired-end alignments
 355     string name = a.Name;
 356     if ( a.IsPaired() )
 357         name.append( (a.IsFirstMate() ? "/1" : "/2") );
 358
 359     // handle reverse strand alignment - bases & qualities
 360     string qualities = a.Qualities;
 361     string sequence  = a.QueryBases;
 362     if ( a.IsReverseStrand() ) {
 363         Utilities::Reverse(qualities);
 364         Utilities::ReverseComplement(sequence);
 365     }
 366
 367     // write to output stream
 368     m_out << "@" << name << endl
 369           << sequence    << endl
 370           << "+"         << endl
 371           << qualities   << endl;
 372 }
 373
 374 // print BamAlignment in JSON format
 375 void ConvertTool::ConvertToolPrivate::PrintJson(const BamAlignment& a) {
 376
 377     // write name & alignment flag
 378     m_out << "{\"name\":\"" << a.Name << "\",\"alignmentFlag\":\"" << a.AlignmentFlag << "\",";
 379
 380     // write reference name
 381     if ( (a.RefID >= 0) && (a.RefID < (int)m_references.size()) )
 382         m_out << "\"reference\":\"" << m_references[a.RefID].RefName << "\",";
 383
 384     // write position & map quality
 385     m_out << "\"position\":" << a.Position+1 << ",\"mapQuality\":" << a.MapQuality << ",";
 386
 387     // write CIGAR
 388     const vector<CigarOp>& cigarData = a.CigarData;
 389     if ( !cigarData.empty() ) {
 390         m_out << "\"cigar\":[";
 391         vector<CigarOp>::const_iterator cigarBegin = cigarData.begin();
 392         vector<CigarOp>::const_iterator cigarIter  = cigarBegin;
 393         vector<CigarOp>::const_iterator cigarEnd   = cigarData.end();
 394         for ( ; cigarIter != cigarEnd; ++cigarIter ) {
 395             const CigarOp& op = (*cigarIter);
 396             if (cigarIter != cigarBegin)
 397                 m_out << ",";
 398             m_out << "\"" << op.Length << op.Type << "\"";
 399         }
 400         m_out << "],";
 401     }
 402
 403     // write mate reference name, mate position, & insert size
 404     if ( a.IsPaired() && (a.MateRefID >= 0) && (a.MateRefID < (int)m_references.size()) ) {
 405         m_out << "\"mate\":{"
 406               << "\"reference\":\"" << m_references[a.MateRefID].RefName << "\","
 407               << "\"position\":" << a.MatePosition+1
 408               << ",\"insertSize\":" << a.InsertSize << "},";
 409     }
 410
 411     // write sequence
 412     if ( !a.QueryBases.empty() )
 413         m_out << "\"queryBases\":\"" << a.QueryBases << "\",";
 414
 415     // write qualities
 416     if ( !a.Qualities.empty() && a.Qualities.at(0) != (char)0xFF ) {
 417         string::const_iterator s = a.Qualities.begin();
 418         m_out << "\"qualities\":[" << static_cast<short>(*s) - 33;
 419         ++s;
 420         for ( ; s != a.Qualities.end(); ++s )
 421             m_out << "," << static_cast<short>(*s) - 33;
 422         m_out << "],";
 423     }
 424
 425     // write alignment's source BAM file
 426     m_out << "\"filename\":\"" << a.Filename << "\",";
 427
 428     // write tag data
 429     const char* tagData = a.TagData.c_str();
 430     const size_t tagDataLength = a.TagData.length();
 431     size_t index = 0;
 432     if ( index < tagDataLength ) {
 433
 434         m_out << "\"tags\":{";
 435
 436         while ( index < tagDataLength ) {
 437
 438             if ( index > 0 )
 439                 m_out << ",";
 440
 441             // write tag name
 442             m_out << "\"" << a.TagData.substr(index, 2) << "\":";
 443             index += 2;
 444
 445             // get data type
 446             char type = a.TagData.at(index);
 447             ++index;
 448             switch ( type ) {
 449                 case (Constants::BAM_TAG_TYPE_ASCII) :
 450                     m_out << "\"" << tagData[index] << "\"";
 451                     ++index;
 452                     break;
 453
 454                 case (Constants::BAM_TAG_TYPE_INT8) :
 455                     // force value into integer-type (instead of char value)
 456                     m_out << static_cast<int16_t>(tagData[index]);
 457                     ++index;
 458                     break;
 459
 460                 case (Constants::BAM_TAG_TYPE_UINT8) :
 461                     // force value into integer-type (instead of char value)
 462                     m_out << static_cast<uint16_t>(tagData[index]);
 463                     ++index;
 464                     break;
 465
 466                 case (Constants::BAM_TAG_TYPE_INT16) :
 467                     m_out << BamTools::UnpackSignedShort(&tagData[index]);
 468                     index += sizeof(int16_t);
 469                     break;
 470
 471                 case (Constants::BAM_TAG_TYPE_UINT16) :
 472                     m_out << BamTools::UnpackUnsignedShort(&tagData[index]);
 473                     index += sizeof(uint16_t);
 474                     break;
 475
 476                 case (Constants::BAM_TAG_TYPE_INT32) :
 477                     m_out << BamTools::UnpackSignedInt(&tagData[index]);
 478                     index += sizeof(int32_t);
 479                     break;
 480
 481                 case (Constants::BAM_TAG_TYPE_UINT32) :
 482                     m_out << BamTools::UnpackUnsignedInt(&tagData[index]);
 483                     index += sizeof(uint32_t);
 484                     break;
 485
 486                 case (Constants::BAM_TAG_TYPE_FLOAT) :
 487                     m_out << BamTools::UnpackFloat(&tagData[index]);
 488                     index += sizeof(float);
 489                     break;
 490
 491                 case (Constants::BAM_TAG_TYPE_HEX)    :
 492                 case (Constants::BAM_TAG_TYPE_STRING) :
 493                     m_out << "\"";
 494                     while (tagData[index]) {
 495                         if (tagData[index] == '\"')
 496                             m_out << "\\\""; // escape for json
 497                         else
 498                             m_out << tagData[index];
 499                         ++index;
 500                     }
 501                     m_out << "\"";
 502                     ++index;
 503                     break;
 504             }
 505
 506             if ( tagData[index] == '\0')
 507                 break;
 508         }
 509
 510         m_out << "}";
 511     }
 512
 513     m_out << "}" << endl;
 514 }
 515
 516 // print BamAlignment in SAM format
 517 void ConvertTool::ConvertToolPrivate::PrintSam(const BamAlignment& a) {
 518
 519     // tab-delimited
 520     // <QNAME> <FLAG> <RNAME> <POS> <MAPQ> <CIGAR> <MRNM> <MPOS> <ISIZE> <SEQ> <QUAL> [ <TAG>:<VTYPE>:<VALUE> [...] ]
 521
 522     // write name & alignment flag
 523    m_out << a.Name << "\t" << a.AlignmentFlag << "\t";
 524
 525     // write reference name
 526     if ( (a.RefID >= 0) && (a.RefID < (int)m_references.size()) )
 527         m_out << m_references[a.RefID].RefName << "\t";
 528     else
 529         m_out << "*\t";
 530
 531     // write position & map quality
 532     m_out << a.Position+1 << "\t" << a.MapQuality << "\t";
 533
 534     // write CIGAR
 535     const vector<CigarOp>& cigarData = a.CigarData;
 536     if ( cigarData.empty() ) m_out << "*\t";
 537     else {
 538         vector<CigarOp>::const_iterator cigarIter = cigarData.begin();
 539         vector<CigarOp>::const_iterator cigarEnd  = cigarData.end();
 540         for ( ; cigarIter != cigarEnd; ++cigarIter ) {
 541             const CigarOp& op = (*cigarIter);
 542             m_out << op.Length << op.Type;
 543         }
 544         m_out << "\t";
 545     }
 546
 547     // write mate reference name, mate position, & insert size
 548     if ( a.IsPaired() && (a.MateRefID >= 0) && (a.MateRefID < (int)m_references.size()) ) {
 549         if ( a.MateRefID == a.RefID )
 550             m_out << "=\t";
 551         else
 552            m_out << m_references[a.MateRefID].RefName << "\t";
 553         m_out << a.MatePosition+1 << "\t" << a.InsertSize << "\t";
 554     }
 555     else
 556         m_out << "*\t0\t0\t";
 557
 558     // write sequence
 559     if ( a.QueryBases.empty() )
 560         m_out << "*\t";
 561     else
 562         m_out << a.QueryBases << "\t";
 563
 564     // write qualities
 565     if ( a.Qualities.empty() || (a.Qualities.at(0) == (char)0xFF) )
 566         m_out << "*";
 567     else
 568         m_out << a.Qualities;
 569
 570     // write tag data
 571     const char* tagData = a.TagData.c_str();
 572     const size_t tagDataLength = a.TagData.length();
 573
 574     size_t index = 0;
 575     while ( index < tagDataLength ) {
 576
 577         // write tag name
 578         string tagName = a.TagData.substr(index, 2);
 579         m_out << "\t" << tagName << ":";
 580         index += 2;
 581
 582         // get data type
 583         char type = a.TagData.at(index);
 584         ++index;
 585         switch ( type ) {
 586             case (Constants::BAM_TAG_TYPE_ASCII) :
 587                 m_out << "A:" << tagData[index];
 588                 ++index;
 589                 break;
 590
 591             case (Constants::BAM_TAG_TYPE_INT8) :
 592                 // force value into integer-type (instead of char value)
 593                 m_out << "i:" << static_cast<int16_t>(tagData[index]);
 594                 ++index;
 595                 break;
 596
 597             case (Constants::BAM_TAG_TYPE_UINT8) :
 598                 // force value into integer-type (instead of char value)
 599                 m_out << "i:" << static_cast<uint16_t>(tagData[index]);
 600                 ++index;
 601                 break;
 602
 603             case (Constants::BAM_TAG_TYPE_INT16) :
 604                 m_out << "i:" << BamTools::UnpackSignedShort(&tagData[index]);
 605                 index += sizeof(int16_t);
 606                 break;
 607
 608             case (Constants::BAM_TAG_TYPE_UINT16) :
 609                 m_out << "i:" << BamTools::UnpackUnsignedShort(&tagData[index]);
 610                 index += sizeof(uint16_t);
 611                 break;
 612
 613             case (Constants::BAM_TAG_TYPE_INT32) :
 614                 m_out << "i:" << BamTools::UnpackSignedInt(&tagData[index]);
 615                 index += sizeof(int32_t);
 616                 break;
 617
 618             case (Constants::BAM_TAG_TYPE_UINT32) :
 619                 m_out << "i:" << BamTools::UnpackUnsignedInt(&tagData[index]);
 620                 index += sizeof(uint32_t);
 621                 break;
 622
 623             case (Constants::BAM_TAG_TYPE_FLOAT) :
 624                 m_out << "f:" << BamTools::UnpackFloat(&tagData[index]);
 625                 index += sizeof(float);
 626                 break;
 627
 628             case (Constants::BAM_TAG_TYPE_HEX)    : // fall-through
 629             case (Constants::BAM_TAG_TYPE_STRING) :
 630                 m_out << type << ":";
 631                 while (tagData[index]) {
 632                     m_out << tagData[index];
 633                     ++index;
 634                 }
 635                 ++index;
 636                 break;
 637         }
 638
 639         if ( tagData[index] == '\0' )
 640             break;
 641     }
 642
 643     m_out << endl;
 644 }
 645
 646 // Print BamAlignment in YAML format
 647 void ConvertTool::ConvertToolPrivate::PrintYaml(const BamAlignment& a) {
 648
 649     // write alignment name
 650     m_out << "---" << endl;
 651     m_out << a.Name << ":" << endl;
 652
 653     // write alignment data
 654     m_out << "   " << "AlndBases: "     << a.AlignedBases << endl;
 655     m_out << "   " << "Qualities: "     << a.Qualities << endl;
 656     m_out << "   " << "Name: "          << a.Name << endl;
 657     m_out << "   " << "Length: "        << a.Length << endl;
 658     m_out << "   " << "TagData: "       << a.TagData << endl;
 659     m_out << "   " << "RefID: "         << a.RefID << endl;
 660     m_out << "   " << "RefName: "       << m_references[a.RefID].RefName << endl;
 661     m_out << "   " << "Position: "      << a.Position << endl;
 662     m_out << "   " << "Bin: "           << a.Bin << endl;
 663     m_out << "   " << "MapQuality: "    << a.MapQuality << endl;
 664     m_out << "   " << "AlignmentFlag: " << a.AlignmentFlag << endl;
 665     m_out << "   " << "MateRefID: "     << a.MateRefID << endl;
 666     m_out << "   " << "MatePosition: "  << a.MatePosition << endl;
 667     m_out << "   " << "InsertSize: "    << a.InsertSize << endl;
 668     m_out << "   " << "Filename: "      << a.Filename << endl;
 669
 670     // write Cigar data
 671     const vector<CigarOp>& cigarData = a.CigarData;
 672     if ( !cigarData.empty() ) {
 673         m_out << "   " <<  "Cigar: ";
 674         vector<CigarOp>::const_iterator cigarBegin = cigarData.begin();
 675         vector<CigarOp>::const_iterator cigarIter  = cigarBegin;
 676         vector<CigarOp>::const_iterator cigarEnd   = cigarData.end();
 677         for ( ; cigarIter != cigarEnd; ++cigarIter ) {
 678             const CigarOp& op = (*cigarIter);
 679             m_out << op.Length << op.Type;
 680         }
 681         m_out << endl;
 682     }
 683 }
 684
 685 bool ConvertTool::ConvertToolPrivate::RunPileupConversion(BamMultiReader* reader) {
 686
 687     // check for valid BamMultiReader
 688     if ( reader == 0 ) return false;
 689
 690     // set up our pileup format 'visitor'
 691     ConvertPileupFormatVisitor* v = new ConvertPileupFormatVisitor(m_references,
 692                                                                    m_settings->FastaFilename,
 693                                                                    m_settings->IsPrintingPileupMapQualities,
 694                                                                    &m_out);
 695
 696     // set up PileupEngine
 697     PileupEngine pileup;
 698     pileup.AddVisitor(v);
 699
 700     // iterate through data
 701     BamAlignment al;
 702     while ( reader->GetNextAlignment(al) )
 703         pileup.AddAlignment(al);
 704     pileup.Flush();
 705
 706     // clean up
 707     delete v;
 708     v = 0;
 709
 710     // return success
 711     return true;
 712 }
 713
 714 // ---------------------------------------------
 715 // ConvertTool implementation
 716
 717 ConvertTool::ConvertTool(void)
 718     : AbstractTool()
 719     , m_settings(new ConvertSettings)
 720     , m_impl(0)
 721 {
 722     // set program details
 723     Options::SetProgramInfo("bamtools convert", "converts BAM to a number of other formats",
 724                             "-format <FORMAT> [-in <filename> -in <filename> ... | -list <filelist>] [-out <filename>] [-region <REGION>] [format-specific options]");
 725
 726     // set up options
 727     OptionGroup* IO_Opts = Options::CreateOptionGroup("Input & Output");
 728     Options::AddValueOption("-in",     "BAM filename", "the input BAM file(s)", "", m_settings->HasInput,   m_settings->InputFiles,     IO_Opts, Options::StandardIn());
 729     Options::AddValueOption("-list",   "filename", "the input BAM file list, one line per file", "", m_settings->HasInputFilelist,  m_settings->InputFilelist, IO_Opts);
 730     Options::AddValueOption("-out",    "BAM filename", "the output BAM file",   "", m_settings->HasOutput,  m_settings->OutputFilename, IO_Opts, Options::StandardOut());
 731     Options::AddValueOption("-format", "FORMAT", "the output file format - see README for recognized formats", "", m_settings->HasFormat, m_settings->Format, IO_Opts);
 732     Options::AddValueOption("-region", "REGION", "genomic region. Index file is recommended for better performance, and is used automatically if it exists. See \'bamtools help index\' for more details on creating one", "", m_settings->HasRegion, m_settings->Region, IO_Opts);
 733
 734     OptionGroup* PileupOpts = Options::CreateOptionGroup("Pileup Options");
 735     Options::AddValueOption("-fasta", "FASTA filename", "FASTA reference file", "", m_settings->HasFastaFilename, m_settings->FastaFilename, PileupOpts);
 736     Options::AddOption("-mapqual", "print the mapping qualities", m_settings->IsPrintingPileupMapQualities, PileupOpts);
 737
 738     OptionGroup* SamOpts = Options::CreateOptionGroup("SAM Options");
 739     Options::AddOption("-noheader", "omit the SAM header from output", m_settings->IsOmittingSamHeader, SamOpts);
 740 }
 741
 742 ConvertTool::~ConvertTool(void) {
 743
 744     delete m_settings;
 745     m_settings = 0;
 746
 747     delete m_impl;
 748     m_impl = 0;
 749 }
 750
 751 int ConvertTool::Help(void) {
 752     Options::DisplayHelp();
 753     return 0;
 754 }
 755
 756 int ConvertTool::Run(int argc, char* argv[]) {
 757
 758     // parse command line arguments
 759     Options::Parse(argc, argv, 1);
 760
 761     // initialize ConvertTool with settings
 762     m_impl = new ConvertToolPrivate(m_settings);
 763
 764     // run ConvertTool, return success/fail
 765     if ( m_impl->Run() )
 766         return 0;
 767     else
 768         return 1;
 769 }
 770
 771 // ---------------------------------------------
 772 // ConvertPileupFormatVisitor implementation
 773
 774 ConvertPileupFormatVisitor::ConvertPileupFormatVisitor(const RefVector& references,
 775                                                        const string& fastaFilename,
 776                                                        const bool isPrintingMapQualities,
 777                                                        ostream* out)
 778     : PileupVisitor()
 779     , m_hasFasta(false)
 780     , m_isPrintingMapQualities(isPrintingMapQualities)
 781     , m_out(out)
 782     , m_references(references)
 783 {
 784     // set up Fasta reader if file is provided
 785     if ( !fastaFilename.empty() ) {
 786
 787         // check for FASTA index
 788         string indexFilename = "";
 789         if ( Utilities::FileExists(fastaFilename + ".fai") )
 790             indexFilename = fastaFilename + ".fai";
 791
 792         // open FASTA file
 793         if ( m_fasta.Open(fastaFilename, indexFilename) )
 794             m_hasFasta = true;
 795     }
 796 }
 797
 798 ConvertPileupFormatVisitor::~ConvertPileupFormatVisitor(void) {
 799     // be sure to close Fasta reader
 800     if ( m_hasFasta ) {
 801         m_fasta.Close();
 802         m_hasFasta = false;
 803     }
 804 }
 805
 806 void ConvertPileupFormatVisitor::Visit(const PileupPosition& pileupData ) {
 807
 808     // skip if no alignments at this position
 809     if ( pileupData.PileupAlignments.empty() ) return;
 810
 811     // retrieve reference name
 812     const string& referenceName = m_references[pileupData.RefId].RefName;
 813     const int& position   = pileupData.Position;
 814
 815     // retrieve reference base from FASTA file, if one provided; otherwise default to 'N'
 816     char referenceBase('N');
 817     if ( m_hasFasta && (pileupData.Position < m_references[pileupData.RefId].RefLength) ) {
 818         if ( !m_fasta.GetBase(pileupData.RefId, pileupData.Position, referenceBase ) ) {
 819             cerr << "bamtools convert ERROR: pileup conversion - could not read reference base from FASTA file" << endl;
 820             return;
 821         }
 822     }
 823
 824     // get count of alleles at this position
 825     const int numberAlleles = pileupData.PileupAlignments.size();
 826
 827     // -----------------------------------------------------------
 828     // build strings based on alleles at this positionInAlignment
 829
 830     stringstream bases("");
 831     stringstream baseQualities("");
 832     stringstream mapQualities("");
 833
 834     // iterate over alignments at this pileup position
 835     vector<PileupAlignment>::const_iterator pileupIter = pileupData.PileupAlignments.begin();
 836     vector<PileupAlignment>::const_iterator pileupEnd  = pileupData.PileupAlignments.end();
 837     for ( ; pileupIter != pileupEnd; ++pileupIter ) {
 838         const PileupAlignment pa = (*pileupIter);
 839         const BamAlignment& ba = pa.Alignment;
 840
 841         // if beginning of read segment
 842         if ( pa.IsSegmentBegin )
 843             bases << '^' << ( ((int)ba.MapQuality > 93) ? (char)126 : (char)((int)ba.MapQuality+33) );
 844
 845         // if current base is not a DELETION
 846         if ( !pa.IsCurrentDeletion ) {
 847
 848             // get base at current position
 849             char base = ba.QueryBases.at(pa.PositionInAlignment);
 850
 851             // if base matches reference
 852             if ( base == '=' ||
 853                  toupper(base) == toupper(referenceBase) ||
 854                  tolower(base) == tolower(referenceBase) )
 855             {
 856                 base = ( ba.IsReverseStrand() ? ',' : '.' );
 857             }
 858
 859             // mismatches reference
 860             else base = ( ba.IsReverseStrand() ? tolower(base) : toupper(base) );
 861
 862             // store base
 863             bases << base;
 864
 865             // if next position contains insertion
 866             if ( pa.IsNextInsertion ) {
 867                 bases << '+' << pa.InsertionLength;
 868                 for (int i = 1; i <= pa.InsertionLength; ++i) {
 869                     char insertedBase = (char)ba.QueryBases.at(pa.PositionInAlignment + i);
 870                     bases << (ba.IsReverseStrand() ? (char)tolower(insertedBase) : (char)toupper(insertedBase) );
 871                 }
 872             }
 873
 874             // if next position contains DELETION
 875             else if ( pa.IsNextDeletion ) {
 876                 bases << '-' << pa.DeletionLength;
 877                 for (int i = 1; i <= pa.DeletionLength; ++i) {
 878                     char deletedBase('N');
 879                     if ( m_hasFasta && (pileupData.Position+i < m_references[pileupData.RefId].RefLength) ) {
 880                         if ( !m_fasta.GetBase(pileupData.RefId, pileupData.Position+i, deletedBase ) ) {
 881                             cerr << "bamtools convert ERROR: pileup conversion - could not read reference base from FASTA file" << endl;
 882                             return;
 883                         }
 884                     }
 885                     bases << (ba.IsReverseStrand() ? (char)tolower(deletedBase) : (char)toupper(deletedBase) );
 886                 }
 887             }
 888         }
 889
 890         // otherwise, DELETION
 891         else bases << '*';
 892
 893         // if end of read segment
 894         if ( pa.IsSegmentEnd )
 895             bases << '$';
 896
 897         // store current base quality
 898         baseQualities << ba.Qualities.at(pa.PositionInAlignment);
 899
 900         // save alignment map quality if desired
 901         if ( m_isPrintingMapQualities )
 902             mapQualities << ( ((int)ba.MapQuality > 93) ? (char)126 : (char)((int)ba.MapQuality+33) );
 903     }
 904
 905     // ----------------------
 906     // print results
 907
 908     // tab-delimited
 909     // <refName> <1-based pos> <refBase> <numberAlleles> <bases> <qualities> [mapQuals]
 910
 911     const string TAB = "\t";
 912     *m_out << referenceName       << TAB
 913            << position + 1        << TAB
 914            << referenceBase       << TAB
 915            << numberAlleles       << TAB
 916            << bases.str()         << TAB
 917            << baseQualities.str() << TAB
 918            << mapQualities.str()  << endl;
 919 }