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