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